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1.0 SUMMARY 


The development of an advanced supersonic panel method is described in detail. The 
basic integral equations of linearized supersonic theory are derived with a discussion on 
boundary conditions providing uniqueness of the solution. Because of the success of the 
method of Johnson and Rubbert for subsonic flow, their geometry and spline system was 
first utilized by essentially replacing the subsonic aerodynamic influence coefficients by 
the equivalent supersonic relations. 

The source method first was tested on a variety of regular three dimensional bodies, 
including circular and elliptic cones and pointed bodies of revolution. The calculated results 
for circular cones and bodies of revolution were in agreement with exact theory; and for ellip- 
tic cones, good agreement with second order linear theory was found. The 6 parameter doub- 
let panel method was tested on cambered wings with zero thickness and linearized boundary 
conditions and gave excellent agreement with exact linearized solutions. 

Using source paneling on both surfaces of thin wings and on inlet nacelles lead to the 
discovery that internal wave reflections induced severe perturbations on the exterior pressure 
distributions. These perturbations were eliminated by using a combined source and doublet 
panel system with interior boundary conditions specified to eliminate the internal perturba- 
tion flow. 

An analysis of the aerodynamic influence coefficients indicated that discontinuities of 
the doublet strength, or of geometry at panel edges, introduce infinite square root singulari- 
ties on the Mach cones emanating from panel corners. For certain panel configurations these 
singularities produced large oscillations in the pressure, even though the doublet strength was 
very nearly continuous across panel edges. On the first attempt to eliminate tills problem, 
the line integrals along panel edges which involve the doublet strength and which contain 
the strongest singularities were discarded; but the method proved inaccurate, indicating a 
strong sensitivity to continuity of doublet strength. 

It became apparent that to insure sufficient accuracy and stability the supersonic panel 
method must maintain continuity of both the doublet strength and geometry. A panel sys- 
tem with all contiguous edges was obtained by dividing the basic four point non-planar 
panel into eight triangular subpanels. A quadratic doublet distribution is applied over each 
triangular subpanel in such a way that the doublet strength is continuous at panel edges, 
leading to a new nine parameter spline for the complete panel. This new spline, with combined 
source and doublet panels to eliminate internal perturbations, was used to compute the flow 
over three wing-body models and yielded pressure coefficients in good agreement with experi- 
ment. The present method appears to be insensitive to how the configuration is paneled, pro- 
vided no gaps occur at panel edges and no subpanels representing solid surfaces are inclined 
at angles with respect to the free stream greater than the Mach angle. 

To close inlets on nacelles for the purpose of eliminating internal flows, superinclined 
panels were developed and tested on a simple nacelle. The superinclined network was capable 
of absorbing the internal perturbations from the lip of the nacelle. With the addition of the 
superinclined networks, the panel method is now capable of computing the flow over wing 
and body combinations which include engine inlets as well. 

For those readers interested only in a simple overall view of the method and the results, 
section 2 the introduction, section 3 on theory, section 4 on the description of the panel 
method and sections 1 2 and 1 3 on results and conclusions may be understood without the 
details provided by the remaining sections. 



2.0 INTRODUCTION 


2.1 HISTORICAL DEVELOPMENT OF PANEL METHODS 

Methods based on linearized theory with singularity distributions over panels represen- 
ting aircraft surfaces have been found especially useful tools for analyzing the aerodynamic 
forces on aircraft. These methods can treat configurations of general shape which are not 
tractible by direct mathematical analysis. Boundary conditions are applied at discrete points 
associated with each panel of the surface. The required integrals are usually evaluated in 
closed form, and a set of linear equations results to be solved for the required parameters. 
These methods have been especially successful for subsonic or incompressible flow and since 
a good summary is given in Ashley and Rodden [ 1 ] and Rubbert and Saaris [2] , they will 
not be discussed here. The most recent method is that of Johnson and Rubbert [3] . Much 
of the technique for supersonic flow presented here is derived from their method. 

In recent years, there are three panel methods based on linearized supersonic flow 
which are noteworthy. Woodward [4] used constant pressure and constant source panels 
and applied tangential flow boundary conditions on wings and bodies. His method was later 
improved by utilizing constant line sources and vortices which vary linearly in the streamwise 
direction on each panel (reference [5] ). 

To obtain differentiability of the computed pressure, Mercer, Weber, and Lesferd [6] 
used singularity splines providing continuity of the vorticity. Difficulties arose at discon- 
tinuities in the plan form for the vorticity method which were alleviated by rounding such 
corners. Piecewise linear pressure distributions were obtained for planar wings with linear 
boundary conditions instead of the stepwise distribution from the constant singularity panel 
methods. 

Morino, Chen, and Suciu [7] describe a method using Green’s formula for the potential. 
The boundary conditions of tangential flow are not applied directly to the potential in the 
useful way but are inserted into the source term in Green’s theorem. Requiring continuity 
of tangential flow near the surface insures that the boundary conditions are satisfied. The 
configuration is divided into quadrilaterals defined by an array of grid points on the surface. 
The surface is approximated for each panel by fitting a hyperboloidal surface through the 
four points which maintains geometric continuity with adjacent panels. The quantities, such 
as values of doublet and source strengths, were assumed constant on each panel leading to an 
influence coefficient method with a set of linear equations to be solved for the values of the 
unknown doublet strengths. 

In the appendix to reference [6] , Mercer, Weber, and Lesferd suggest a doublet Mach 
line panel method using a fourth degree polynomial in characteristic coordinates for the 
doublet distribution on each panel to analyze planar wings. This type of paneling has a num- 
ber of positive features: 

1 . Discontinuities in pressure occur across Mach lines from planform corners and can be 

easily taken into account with Special Mach line paneling. 
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The aerodynamic influence coefficients in the characteristic coordinates are very sim- 
ple, closed form expressions. 

3. Because of the domain of influence of the characteristic strips, the matrix of coefficients 
for the parameters can be made triangular by proper ordering, and hence the solution is 
exceedingly fast. 

4. A continuous pressure distribution is obtained. 

This method, along with a similar source panel technique, was derived and tested as part 
of the present contract. The basic theory and results are presented in a separate contractors 
report (reference [8] ). The source panel method was stable and accurate for both analysis 
and design boundary conditions. The doublet method was also stable for the region down- 
stream of the Mach line from the corner formed by the supersonic leading edge but was unstable 
for the region downstream of the Mach line from the comer formed by the supersonic and 
subsonic portions of the leading edge. The cause of this instability was not pursued, principally 
because the method was confined to planar wings with linearized boundary conditions and 
hence has limited application, but also because the subsonic method of Johnson and Rubbert 
[3] appeared to be a more promising approach. 

The most advanced panel method for subsonic flow is that developed by Johnson and 
Rubbert [3] . The technique uses a quadratic approximation for curved panels with a quad- 
ratic distribution of doublet strength and a linear distribution of source strength applied to 
each panel. The vanishing of the normal component of the mass flux is applied as boundary 
conditions to solid surfaces instead of the vanishing of the normal component of velocity as 
in most of the earlier methods. This is discussed in section 3. The aerodynamic influence 
coefficients are correct to the first order in relative panel curvature. The method has the 
following properties: 

1 . It is insensitive to how the configuration is paneled; and, hence, allows the use of auto- 
matic paneling programs. It does not require special experience in applying the program 
to practical engineering problems. 

2. It is economical in computation cost. For a given accuracy, the method requires fewer 
panels than the constant singularity panel techniques. 

3. It offers a wide variety in the choice of modeling techniques, including both design and 
analysis boundary conditions. 

The numerical method is stable, accurate, flexible, and efficient, essential properties 
for any numerical method to be useful and practical for engineering design and analysis. It 
was felt that the approach was ideal for supersonic flow as well. The techniques for treating 
the geometry and the mode of paneling in the subsonic method were utilized with very little 
change in the initial supersonic method, along with the quadratic doublet and linear source 
distribution on each panel. In subsonic flow, disturbances due to discontinuities of the doub- 
let strength at panel edges or gaps in panel geometry decrease in intensity with distance from 
their origin. In supersonic flow, similar discontinuities produce disturbances which do not 
decay with distance but introduce infinite singularities along Mach cones emanating from panel 
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comers. Hence, continuity requirements for supersonic doublet panels are considerably more 
strict than for subsonic. To overcome this difficulty, a paneling system with contiguous 
edges and a doublet spline with continuous doublet strength were developed. 


2.2 TECHNICAL APPROACH AND DEVELOPMENT OF THE PRESENT 
SUPERSONIC PANEL METHOD 


With the considerable success of the subsonic panel method of Johnson and Rubbert 
[3] , it was natural that the same method of paneling and the same splines be utilized directly 
with the supersonic aerodynamic influence coefficients replacing the subsonic ones. The 
basic theory for linearized supersonic flow is presented in section 3. The perturbation poten- 
tial is derived as an integral over the surface of sources and doublets in section 3.7 in the form 


0(R) = - 


27T 



o(R n ) ds 

R — + 
R B 


b2 

27 r 



(R 0 ) (Rq-R) . nds 
R B 3 


where the asterisk indicates finite part (see section 3.4), o is the source strength, /i the doublet 
strength, 

r b = V (x-xq,) 2 - B 2 (y-yo) 2 - B 2 (z-z Q ) 2 


is the hyperbolic distance, R = (x,y,z), Rq, the integration point, and B 
= (n x , n y , n z ) is the unit normal vector to the surface. 


v 




and n 


The boundary conditions that the normal component of the linearized mass flux vector 
vanish on solid surfaces is given by the combination of equations (3.16), (3.12), and (3.8) and 
takes the form 


(1-B^0 X ) n x + 0y n y + 0z n z ~ 0 


for the freestream velocity in the x direction. This is the basic boundary condition applied 
to surfaces over which the supersonic flow is to be determined. 

The configuration to be analyzed is divided into quadrilateral panels by a set of grid 
points on the surface. As in reference [3] , an average plane is defined by the four midpoints 
of the lines joining the grid points. A linear source and a quadratic doublet distribution are 
prescribed over the projection of the panel onto the average plane and the integrals were 
evaluated in closed form. Associated with each panel, are a value of the source and a value 
of the doublet strength together with two appropriate boundary conditions. This leads to a 
set of linear equations to be solved in the form 
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n 

ayXj — bj , j 1,2,. . n 

i= 1 

where ay is the matrix of aerodynamic influence coefficients, Xj are the singularity strengths. 
Once the Xj are found the velocity and pressure may be computed anwhere in the flow field. 

Johnson and Rubbert [3] also included a correction to the aerodynamic influence coef- 
ficients due to panel curvature, by expanding the potential for a curved panel in terms of the 
small quantities defining the panel shape and retaining only the first order terms. The same 
approach failed for supersonic flow because of the singularity on the Mach cone which was 
of higher order for the curvature terms and, hence, the pressures were non-integrable. With 
flat panels, the method produced excellent results for a wide range of configurations, and 
these are reported in section 12. 

Hess [9] has shown that a discontinuity in doublet strength induces the same velocity 
field as a line vortex. (The relationship between doublet sheets, vortex sheets and line vor- 
tices is explained in Appendix A of reference [9] ). Even when the doublet is continuous from 
one panel edge to the adjoining panel edge, gaps between panels will produce infinite singu- 
larities on the Mach cones from the panel corners because the gap prevents the cancelling of 
the two line vortices. These singularities in supersonic flow do not decay with distance from 
the edge as in subsonic flow but propagate along Mach cones nearly unattenuated (see section 
A8.1 in Appendix A). For some paneling configurations, these disturbances impinge upon 
control points and cause large oscillations in pressure. It was felt that since the discontinui- 
ties in doublet strength were small, the problem could be eliminated by discarding the line 
integrals of the doublet strength along panel edges which produce the line vortices. (For 
example, the second integral of the equation following equation (7.1)). When this was done 
the method proved to be inaccurate even for small discontinuities of doublet strength or geo- 
metry, although this approach was successful for subsonic flow. 

It became apparent that, for supersonic flow, the gaps at panel edges must be eliminated 
by a panel system that maintains continuity of geometry and by a spline system which ensures 
continuity of doublet strength. This was achieved by dividing the basic non-planar panel 
formed by the four points on the curved surface into 8 triangular panels in such a way that no 
gaps occur in the geometry (see figure 1). A quadratic distribution of doublet strength is pres- 
cribed over each triangular subpanel in such a way that doublet strength is continuous over the 
entire surface. This leads to a 9 parameter spline instead of the former 6 parameter spline; 
and with the combined source and doublet paneling yields satisfactory solutions for the super- 
sonic flow over complicated wing body configurations. 

Section 3 presents a derivation of the basic integral equations of linearized theory with 
appropriate boundary conditions, while section 4 presents a more detailed description of the 
panel method. Derivation of specific formulas for the method are given in Section 5 through 
1 1 . For those interested only in an overall view of the panel method, section 4 and also sec- 
tions 1 2 and 1 3 giving results and conclusion may be read without the necessity of having 
read the complete theory. 
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3.0 THEORY OF LINEARIZED COMPRESSIBLE FLOW 


3.1 DERIVATION OF BASIC EQUATIONS 

To derive the basic equation we follow the analysis of Ward [10] 
equations for steady, inviscid flow in the form 

v x V c = 0 


v . pV c = 0 


pV c • vV c + vp = 0 


and consider Euler’s 


(3.1) 


(3.2) 

(3.3) 


where V c = U + V, U is the free stream velocity vector, and V the perturbation velocity vec- 
tor. The quantities p and p are the pressure and density, respectively. With the subscript 0 
denoting free stream values of the quantities, linearizing of equations (3.1) to (3.3) yields 


vx V = 0 

(3.4) 

O 

< 

• 

<1 

+ 

C| 

• 

< 

*o 

o 

II 

o 

(3.5) 

PqU • vV + v (p-pg) = 0 

(3.6) 


? 2 

Since for isentropic flow dp = c^dp - cq dp, where c is the velocity of sound, eliminating p 
and p between equations (3.5) and (3.6) gives a single differential equation for the velocity, 
namely, 


v.[V-U(U*V)/c 0 2 ] = 0 


or 


V. W = 0 


(3.7) 
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where W = V - U(U • V)/cq 2 and we have made use of the fact that U is a constant vector. 
When we assume that U is a vector in the x direction, W becomes 

W = [(l -M 2 )u,v,w] 

where M 2 = U»U/cq 2 is the square of the free stream Mach number. From equation (3.7) 
we see that the quantity W is conserved in the flow. In the following we shall show that W 
is proportional to the perturbation mass flux vector. 

To find the local linearized pressure we consider Bemouilli’s equation in the form 

•P 


(3.8) 


/ 

Po 


dp/p + v c . V c /2 = U . U/2 


Expanding 1/p in powers of p - pq 

1 /P - 1/po - (dp/dp)o (p-PoVpo 2 + 

= I/p °“fe) ~T~ + -“ 

\ PO / cq z pq 


(3.9) 


and substituting into equation (3.9) yields 


(P-PQ) 

PO 


i /p-pq y 

2c 0 2 V PO ) 


+ . . . = - U . V - V . V/2 


since V c = U + V. This can be solved for (p - Pq)/pq by inverting the series. Thus, retaining 
no terms higher than quadratic gives 


P-PQ /P 0 = -U.V-V.W/2 


(3.10) 


The quantity p Q W is now shown to be the linearized perturbation mass flux vector. Con- 
sider the difference between the local mass flux and that for the free stream. This is 

P V C - POU = (P-PO) u + P0 V 


with 

P-PO = (dp/dp)o (p-Po) + • • • - (p-Po)/ c 0 2 +•• • • 
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we have from equation (3.10) 


pV c -p 0 u = p 0 [V-U(U.V)/c 0 2 + ...] = P 0 W + . . . (3.11) 


Hence equation (3.7) is shown to be the conservation equation for the perturbation mass 
flux vector. 


The conditions of irrotationality, equation (3.1), is satisfied by introducing a velocity poten- 
tial. Let 


V = U 


(3.12) 


then the conservation equation (3.7) becomes, with U in the x direction, 

( 1 -M^) $yy 0zz “ ® 


(3.13) 


and is easily recognized as the classic Prandtl-Glauert equation for linearized compressible 
flow. From equation (3.10) the pressure coefficient becomes 


Cp = 1 — P — =-20 x -(j3 2 0 x 2 + 0y 2 + 0 z 2 ) (3.14) 

^P 0 U 0 2 


where 


P = 1 - M 2 


When j3 2 0 x 2 is dropped this is the slender body approximation to the pressure. For flat wings, 
the quadratic terms are often neglected. However, we often compute pressure by the slender 
body approximation or by the complete isentropic relation. 

3.2 UNIQUENESS OF THE MASS FLUX BOUNDARY CONDITIONS 

IN SUBSONIC FLO V/ 

To demonstrate uniqueness of the solution of equation (3.7) or (3.13) under certain 
boundary conditions for subsonic flow, we apply the divergence theorem to 0W for a volume 
V enclosed by a surface S. We obtain, using v* W = 0, 

f rf>W • n ds = f W.v0dv= f (p^u^ + V 2 + w 2^ dv (3.15) 

4 ■'v *V 
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where n =( n x , n y , n z ) is the outward normal to the surface S. Let S be divided into two 
parts with 0 prescribed on Sj and W . n prescribed on S2. Since the differential equation is 
linear, the difference of two solutions is also a solution and equation (3.15) holds. Let 0 
and W denote the difference of two solutions which satisfy the same prescribed boundary 
conditions. Then the surface integral on the left hand side of the equation (3.15) vanishes. 

The volume integral on the right side is identically equal to zero. Since the terms of the inte- 
grand are always positive, the integral can be zero only if the terms of the integrand are zero 
everywhere, and the two solutions are identical. Hence, the solution of V- W = 0 will be 
unique with 0 or W • n prescribed as boundary conditions on the surface S. This can be shown 
for infinite regions when 0 goes to zero like 1/r as r goes to infinity. In the light of the uni- 
queness theorem it is appropriate in our panel methods for analyzing the flow over bodies to 
use linearized mass flux boundary conditions. On solid boundaries this takes the form 

(U + W) • n = 0 (3.16) 

or for the velocity potential in equation (3.12) and a free stream velocity in the x direction 

( 1 + 0^0 X ) n x + 0 y n y + 0 z n z = 0 (3.17) 

This equation can be written in a different form by introducing the conormal vector 

n c = 0 2 n x , n y , n z 

We obtain for equation (3.1 7); 

n x + 30/3 n c = 0 

The conormal derivative of the potential then can be interpreted as the normal component of 
the perturbation mass flux to the surface. 

The differential equation (3.13), the pressure relation in equation (3.14) and boundary 
conditions in (3.17) comprise the fundamental boundary value problem we are solving for 
analyzing aircraft configurations by the methods described here. When the flow is subsonic, 
the equations can be converted to the incompressible flow by introducing the incompressible 
velocity potential 

0 = 0 i/p (3.18) 
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and scaling the x variable by 


x -> x//3 (3.19) 

Since n x ~* n x //3, we obtain 


W (0ix> $iy/& $iz Iffy 

(3.20) 

^ixx + $iyy + 0izz = 0 

(3.21) 


Cp C pi /^ (3.22) 

Cpi = - 20j x - + 0jy2 + 0j z 

^1 0ix) n x 0iy n y ^iz n z “ ^ (3.23) 


Note that the compressible pressure coefficient including the quadratic terms is related to 
the incompressible value by the factor l/j3^. This was first pointed out by Goethert [11]. 
Equations (3.21) and (3.23) form the classical Neumann boundary value problem for Laplace’s 
equation. Solutions to this problem are unique. This is another proof of uniqueness for the 
mass flux boundary conditions in linearized subsonic flow. Equations (3.21), (3.22), and 
(3.23) are the basic equations for the subsonic panel method described in Johnson and Rubbert 
[3] and which is utilized to a great extent in developing the supersonic method to be presented 
in the following. 

Equation (3.17) and the definition of the conormal hold for supersonic flow as well. 

Since M > 1 , we obtain 


(l -B 2 0 x ) n x + 0y n y + 0 z n z = 0 


where B = s / - 1. The conormal then takes the form 

n c B2n x , ny, n% 
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3.3 EULER MOMENTUM THEOREM FOR MASS FLUX 
BOUNDARY CONDITIONS IN LINEARIZED THEORY 

We shall show that the mass flux boundary conditions are also compatible with Euler’s 
momentum theorem when p, the pressure, is given by equation (3.10). For a closed surface 
S with no internal sources and for which the velocity vector V c is not singular, the theorem 
takes the general form: 


f [pn + p(v c .n) V c ]ds = 0 
For linearized flow the mass flux vector pY c becomes 


pV c - pq (U + W) 


and we have the momentum integral 


/[- + p 0 (U + W) . n (U + V)] ds 


(3.24) 


This can be shown to be valid for any surface S enclosing a volume V of the fluid 
when the pressure is given by equation (3.10). After eliminating the pressure, the momen- 
tum integral in equation (3.24) becomes 


P0 fn . [- I(U • V + V . W/2) + (U + W) (U + V)] ds 


where I is the identity matrix for which n . I = n . Applying the divergence theorem 


f n • A ds = f v • A dv 

•'s 

replaces the surface integral by the volume integral 

) /[-v(U.V)-v(V.W)/2 + U-v W + W. vV + (U + V)v. W]dv 


P 0 
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where we have made use of the fact that U is a constant vector and also 


v. (1^) = v \J/ 


for a scalar and 


V. (v 1 V 2 )=V 2 v . Vj +Vj .V V 2 


for vectors V] and V 2 . 

Now applying the formula from vector analysis 

v, X (v 2 X v 3 ) = (Vj . V 3 ) V 2 - (V! . V 2 ) V 3 

to W X ( v X V) and noting that v operates only on V yields 

Wx(vxV) = v[W] . V - W . v V 


Since the vectors V and W are related by a constant symmetric matrix in the statement 
following equation (3.7) they have the property that 


We finally obtain 


3W 

dx 


V = 


w 


3V 

3x 


Wx (vxV)= v(V • W)/2 - W • vy 
Similarly, as U is a constant vector, we have 


U x (vx V) = v(U • V) - U • W 
With these substitutions, the volume integral finally becomes 


p° 1 J [ - Ux (v x V)- Wx (vx V) + (U. V) v. W] dv 


Since we assume potential flow, then V X V = 0 and V. W = 0. The integrand of the volume 
integral vanishes identically thereby proving that equation (3.24) holds for any arbitrary 
surface enclosing the flow when the pressure is given by equation (3.10). 
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Let S consist of a body S b , both sides of a_doublet sheet S v , and a surface S' enclosing 
the body. It will be shown in section 3.7 that W . ft is continuous across the doublet sheet. 



Cross-section of Surface S 

Since (U + W) • n = 0 on the body and W . n is continuous across the doublet sheet, the 
integral over and S v yields the body force 

F = C pnds 

S b 


Hence the body force F can be found by an integral of the pressure and momentum through 
a surface S' enclosing the body using equation (3.24); namely, 


= - /*[pn + po (U + V) (U + W) . ft] ds 

“o' 


The surface S' is arbitrary but care must be taken in integrating about edges of the doublet 
sheets with singularities in the potential. For greater details the reader may consult Ward 
[10] and the applicable references. 
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3.4 INTEGRAL EQUATIONS FOR SUPERSONIC LINEARIZED FLOW 

Except for the uniqueness proof and scaling of the velocity potential, the foregoing 
analysis is directly applicable to supersonic flow. To derive an integral equation method 
we apply the divergence theorm to the quantity 

0lW 2 -02Wi (3.25) 


where the subscripts 1 and 2 denote two distinct solutions of equation (3.7). Using the 
identity • W 2 = v0 2 • Wj , we obtain 


J\ W 2 - 0 2 Wj) • n ds = J* ^0j v W 2 - 0 2 v Wj )dv (3.26) 


S V 

For the subsonic method 0 2 was chosen as the point source 


<t >2 = i //(*-* o) 2 + 0 2r2 ’ r2 = (y-yo) 2 + ( z_z o) : 


(3.27) 


which satisfies V- W 2 = 0 except at the point xq = x, yQ = y, and zq = z. With 0 2 from 
equation (3.27) in equation (3.26), the surface integral yields a solution for 0j in terms of 
source and doublet distributions over the surface. An analogous procedure was used by 
Hadamard [ 12] , Ward [10], and by Heaslet and Lomax [13] for supersonic flow. Since 
M > 1 , the coefficient of the 0 XX term in equation (3.1 3) is negative and the differential 
equation becomes hyperbolic. The supersonic equivalent to equation (3.27) is then 


02 - 1/R B 


(3.28) 


where 


Rg = (x-xq) ^ _ B2f2 and b2 = M^-1. 

The supersonic case is much more complicated than the subsonic, since the equivalent point 
source has singularities along the surface x - xq = ± Br as well as at the source point xq, 
yQ, zq. With 0 2 substituted into equation (3.26) some of the integrals do not exist in the 
usual mathematical sense and require special treatment first introduced by Hadamard [12] 
and explained briefly farther on in this section of the report. 

It is worthwhile to examine the fundamental differences between subsonic and super- 
sonic flows. When a disturbance occurs at a time t = 0 the wave propagates into the fluid 
at the speed of sound while at the same time the disturbance is being convected by the 
stream. In an infinite time with a subsonic free stream, the disturbance will fill entire 
space (see tig. 2). When the stream is supersonic, the disturbance is convected away from 
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its origin at a speed faster than its propagation rate into the medium. Consequently, the 
disturbance is confined to a conical region formed by the envelope of spherical wave fronts 
of the convected disturbance source as shown in figure 3. The sine of the cone angle is 
the ratio of sound speed eg to the free stream velocity Ug, or the reciprocal of the Mach 
number. The region of influence of the disturbance from a point is easily recognized from 
figure 4 as the region bounded by the downstream cone from ( xg, yg, zg). In this 
region the radicand in the point source 02 in equation (3.28) is positive or 

lx-xg I > Br 


The Mach cone has two nappes. The downstream cone x = xg + Br bounds the domain 
of influence of the point xg, yg, zg, while the upstream cone x = xg - Br bounds the 
region which influences the point xg, yg, zg, or the domain of dependence of the point 

x o> yo> z 0 ( see fi s- 4 )- 

The Mach cone is also a characteristic surface associated with the differential equation. 
Across characteristic surfaces the normal derivative of the perturbation potential may be 
discontinuous while the tangential derivative remains continuous. This can occur only if 
the differential equation can be expressed in terms of derivatives along these surfaces. 
Characteristic surfaces are also the surfaces for which the normal component of the free 
stream velocity is equal to the speed of sound and hence can only exist in supersonic flow. 

These characteristic surfaces can represent shocks in linearized supersonic flow 
since in the limit as shock strength approaches zero, the shock approaches the characteristic 
direction. For exact supersonic flow, entropy increases across a shock and the flow is 
usually no longer irrotational. Since the entropy is of the third order in shock strength, 
the assumption of irrotational flow is still valid to the second order in linearized supersonic 
flow. For flows containing only weak shocks, linearized theory can be expected to yield 
satisfactory results. 

Discontinuous expansion waves also occur at characteristic surfaces in linearized 
theory. In exact theory, the expansion waves are in the form of a continuous centered fan 
of characteristics instead of a single discontinuity. The approximation of expansion waves 
by an ‘expansion shock’ is valid near the surface but is a poor representation of the flow at 
greater distances away from the surface. Linearized theory solutions involving rapid 
expansions can be expected to yield good pressure distributions on the surface but the 
induced flow at large distances from the expansion surface will not be accurately described. 

Because of the inverse square root singularity on the Mach cone of the supersonic 
source, the derivation of the integral equation method for subsonic flow using source and 
doublet surfaces cannot be carried over directly into supersonic flow, since derivatives of 
the source occur. Volterra (see ref. [ 14] ) avoided this difficulty by using a different 
fundamental solution for 0 2 which vanishes on the Mach cone. Hadamard [12] got around 
the difficulty by defining the finite part of the divergent integrals resulting from using 
the acoustic source equivalent to equation (3.28) in solutions of the wave aquation. This 
approach was also applied by Ward [10], Robinson [ 1 5 ] , and Heaslet and Lomax [13]. 

To illustrate this concept, consider the following integral 
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f(x)dx/(x-a)3/2 


If f'(a) exists and is bounded in the interval a to be, then 1(e) exists. To assign a meaning 
for e = 0 we write 1(e) in the form 


/ b y* b 

ff(x) - f(a)] dx/(x-a)3/2 + f(a) J dx/(x-a)3/2 
a+e a+e 


f [f(x) - f(a)] 
a+e 


dx/(x-a)3/2 - 2f(a)/(b-a)l/ 2 + 2f(a)/e 1 / 2 


As e-^0 the first integral exists as an ordinary improper integral. We further write 


I(e) = 



[f(x) - f(a)]/(x-a)3/2 _ 2f(a)/(b-a) 1 /2 + 2f(a)/el/2 



[f(x) - f(a)] dx/(x-a)3/2 


Hadamard [12] defined the first two terms as the finite part of 1(0). Ward [10] used the 
notation *1 to denote this finite part. Hence 

*1 = lim {l(e)-2f(a)/e 1 /2 + 0(e I /2)} 
e~* 0 

The unique finite part of an infinite integral is demonstrated here for a single variable but 
it can be generalized to integration of several variables. The finite part can be shown to 
have the following properties (ref. [ 10] ): 

1 . If I converges then *1 = I. 

2. The value of *1 is invariant under coordinate transformation provided the transformation 
itself is not singular on the surface of singularity of I. 
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3. 


If *1 involves vector quantities, then *1 is invariant under rotation of the coordinate 
axes. 


4. If 1(c) vanishes for all e, then *1 = 0. 


5. Differentiation of finite part integrals with respect to a parameter may be moved 
inside the integral. For example, 



3f/3a ds 


Hadamard [12] also showed that the divergence theorem holds for the finite part 
integrals, i.e., 



. Wdv 



W. nds 


and Robinson [15] extended the principle to include Stokes theorem 



vx V. nds 



V. d£ 


With the concept of finite part of an integral to give meaning to the divergent integrals, 
we are now able to proceed in the same manner as in incompressible flow. We first show 
that the mass flux from <j > 2 in equation (3.28), through a surface enclosing the point xq, 

Yq, zq is independent of the choice of surface. We can assume without loss of generality, 
that xq, yq, z jQ is the origin of coordinates. Let this surface contain the closed cylinder 
-a < X < a, r^ = y2 + < (a/B)2 } which intersects the downstream Mach cone at x = a. 

On this cylinder, 0 2 and W 2 are non-zero only on the plane x = z. 

To evaluate the mass flux we consider the surface S(e) interior to the Mach cone as 
in figure 5. Then the mass flux is given by 

J(e)= J* W • n ds = -B j* 302 /3x dydz = J* B^a dydz/y^a^ - B^r^)^ 
s (<0 s< e ) s< e ) 
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Since r? = a/B - e 


/ a 1 

J(e) = -27r + 27r v / 1 +B 2 e i/ 2 

and 

J* = lim [j(e) - 27r/e 1 / 2 v 4/(l + B2) ] = - 2 tt 
e - * - 0 


This gives the mass flux from the plane cutting the cone. We need to show that for 
two surfaces, S\ and S 2 , which cut the cone for x > 0 as shown in figure 5, the mass flux 
through each surface is the same, or 

W • n ds = ^ W . n ds 

Si S 2 



J\ + B 2 
B 


Applying the divergence theorem to the volume between the two surfaces yields 



W • n ds = 



W . n ds = 



* 


V 


v 


W ds 


where V is the volume enclosed by the two surfaces and the Mach cone. Since this volume 
does not contain the origin where W is not defined, v • W = 0 from equation (3.7) and the 
mass flux is seen to be independent of the choice of surface enclosing the source point. 

The fundamental solution 


0 = - l/27rRg 


(3.29) 


then represents a true supersonic source of unit strength. 

3.5 INTEGRAL EQUATION FOR SUPERINCLINED PANELS 

We now derive an integral equation for superinclined surfaces, i.e., surfaces inclined 
with the free stream at angles greater than the Mach angle. We consider the volume in 
figure 6 enclosed by the surface S, the upstream Mach cone T and the circle Sj, near to but 
excluding the vertex of the cone at the point x, y, z. For 0 2 equation (3.26) we use 
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0 2 = 1/ r B for r B 2 > 0, xq < x 
= 0 at all other points 


and for <j )\ , we simply write 0. Since by equation (3.7)V . W = V . W 2 = 0, equation (3.26) 
yields 


f 

s+r+Sj 


[*w 2 -0 2 w] 


n ds = 0 


The integration over the cone surface T can be eliminated since the integrand is singular 
there and, hence, itjnakes no contribution to the finite part. The integral over the small 
circle Sj involving W 2 becomes 


J' 01 2 • n ds « $ f W 2 . n 
S 1 S 1 


ds 


where the tilda denotes average value over the surface S\ . This is the same integral used 
for computing the mass flux from a point source, but with the normal directed in the 
opposite direction. Hence, we obtain, by letting the surface S| approach the cone vertex, 


^ limit 

Si — 0 


1 


0 W 2 • n ds = 2tt 0(x,y,z) 


(3.30) 


Similarly, for the other integral over S\ 


I 

S] aj u 0 


2tt ~ (x - xq) /B 


rdrdg 


y(x-xo) 2 - BV 


= 2tt 


W . n (x-xq) 


B 


/' 


pdp A 

= 2tt (x-xq) W . ft/B 
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The contribution of Sj from 02 vanishes as x - xq goes to zero. Finally, the integral 
equation for superinclined surfaces takes the form 



f [W. n0 2 -0W 2 . n]ds 


After substituting 0 2 anci 


m> - £ 



W . n ds 

Rb 


2t r 



(3.31) 


It is not necessary to take the finite part of the first integral since the integrand has an 
integrable singularity. Since 


w 2 : 


A 

: n • 


f _ g2 2 _2J 

\ 9xq ’ 9yo’ 9 zQj 


'A (Rq-R) »n 

Rb/ R B 3 


the first integral of equation (3.31) represents a distribution of sources while the second 
integral a distribution of doublets. It is easy to show that 0 takes on the prescribed values 
on the surface S. As the point x, y, z moves to the surface, the integration over S 
becomes the same as the previous integration over near the vertex of the Mach cone, 
but with the normal now directed inwards toward the vertex. Thus, 



02 - 0 W 2 • n) ds 



02 




W 2 • n ds = 0 


Similarly, W • n can be shown to take on the assigned value. 

3.6 UNIQUENESS OF THE SOLUTIONS WITH POTENTIAL AND LINEAR MASS FLUX 
BOUNDARY CONDITIONS ON SUPERINCLINED SURFACES 

It is possible to show that when 0 and W . n are prescribed on a surface inclined at 
an angle greater than the Mach angle, that the solution is unique. To that end, we apply 
the divergence theorem to the matrix quantity 

VW - (V . W) 1/2 
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where I is the identity matrix. We obtain 


J {Vv . W + W. vV- v(V. W)/2} d v =y" ( VW . n - V . Wn/2) ds 


where we have used 


(3.32) 


v (VW) = W v. V + V. vWand v . (V • W) I = v(V • W) 


For irrotational and_source free flow^the vohame integrals disappear since we showed 
in section 3.3 that W. VV - V(V. W)/2 = WxVxV = 0. Consider the x component of 
this surface integral: 


f k ($y n y n z^z ~ ^0x n x) ( ®^x^ + 0y^ ^z^) n x/^J 


ds - 0 


(3.33) 


Simplifying and factoring yields 



(3.34) 


Denote by 0 S the boundary values of 0 and its derivatives. Then the squared terms in the 
first two pairs of parentheses are 0 s ^y and 0 S ^ Z , respectively. Furthermore, if T, say, is a 
characteristic surface, then 


- (n y 2 + n z ^/n x 2 = 0 


On such a characteristic surface T, the integral becomes 

- Si g n ( n x) f [0 S y 2 + ^sz 2 ] d y dz 

r 
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Let S consist of three parts as shown in figure 7. 


1 . The characteristic cone x = xj - B j(y~y 2 + ( z_z l) ^ denoted by I\ 

2. The plane x = 0 with B 2 £(y-y l) 2 + ( z “ z i) 2 J < x j 2 denoted by S \ , and on which 
boundary conditions are to be applied. 

3. The plane x = xq, 0 < xq < xj , with B 2 [(y-yi) 2 + (z-zj) 2 ] < (xq- x i) 2 
denoted by S3. 

Since on Sj and S3, we have n x = -1 and 1, respectively, with n y = n z = 0, then the 
integral in equation (3.34) becomes 

- J f ( <t > sy 2 + hz 2 ) d y dz + If + h 2 + ^ dydz i 

r Sj 

-J J [B20 x 2 + 0 y 2 + 0 y 2 + 0z 2] dydz = 0 (3.35) 

s 3 


since on T, n x > 0. 

Let 0 = 0] - 02 where the subscripts 1 and 2 denote two solutions with the same 
prescribed values of 0 and 0 X on the plane x = 0, i.e., on the surface S] . Then 0 = 0 X = 0 
on S] and the second integral over S\ vanishes in equation (3.35). It follows from the 
form of the remaining integrals that 

^sy - 0sz “ 0x “ ^y ~ 0z “ 0 


The surface S3 may be varied arbitrarily along the cone axis. Hence, 0 X = 0y = 0z = 0 
everywhere inside the cone. Therefore, when 0 X and 0 are prescribed on S \ , then 0 is 
uniquely determined inside the Mach cone from the point X], yj, zj in figure 7. 

Under the section on coordinate transformations, the hyperbolic distance and the 
differential equation are shown to be invariant under a rotation about the x axis and 
under an oblique transformation in the plane containing the x axis. If the transformation 
is chosen for a superinclined panel so that the new variables y, z are in the plane of the 
panel and x is the variable out of the panel; then, for panels inclined to the free stream at 
angles greater than the Mach angle, the area element and the conormal derivative 
transform according to 

(d0/9n c )ds (30/dx)dzdy (3.36) 
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Ward [10] states that for superinclined surfaces, the conditions 


0, 30/3x 


are generalized to 0 and 30/8 n. The boundary conditions should more correctly be 
generalized to 

0, 30/3n c or w • n 


prescribed on the panel. 

3.7 INTEGRAL EQUATION FOR SUBINCLINED PANELS 

To apply panel methods for the solution of -the supersonic flow over wings and 
bodies, we must consider subinclined surfaces, i.e., surfaces inclined at angles less than 
the Mach angle. In section 3.6, it was shown that boundary conditions can only be 
applied on the downwind side of the superinclined surface. For the subinclined surface, 
boundary conditions may be applied to either or both sides. Consider an almost planar sur- 
face. If the surface is slender and lies within the Mach cone, then the perturbed flow is con- 
fined to the downstream Mach cone emanating from its vertex. If the surface cuts the down- 
stream Mach cone emanating from the forwardmost point of the surface, then the region of 
disturbances lies downstream of the characteristic surface formed by the envelope of Mach 
cones emanating from the leading edge (see fig. 8). This is the surface at which the How 
first experiences a disturbance. To apply equation (3.26) we consider 2 volumes, V j and 
V2, illustrated by the cross section drawing in figure 8. The volume V] is bounded by the 
following surfaces: 

Tj consisting of the leading edge characteristic surface cut out by the Mach cone T 
from the point x, y, z. 

S e the circle of radius e formed by the plane excluding the highest order singularity 
at x, y, z from the volume V. 

S w the portion of the subinclined surface cut off by the Mach cone where 
boundary conditions are to be applied. 


We have already shown in equation (3.30) that as e^O the contribution to the surface 
integral S e was equal to 


27r0(x,y,z) 


( 3 . 


The direction of the conormal along a Mach cone or any characteristic surface lies along the 
generator and hence W • n is the perturbation tangential velocity. Because W . n is zero 
ahead of the characteristic surface then W • n is zero on the downstream side since this 
is required by the continuity of mass and momentum across the surface. Hence, 0 is 
constant on the characteristic surface Tj. Since the addition of a constant to the velocity 
potential does not affect the velocity, we may choose this constant as zero to eliminate 


the integral over the surface T j . Since the finite part cancels the integral over the Mach 
cone T where the integrand is infinite, we obtain finally from equation (3.26). 


0(x,y,z) = - T = 


(0 + W 2 - W + 02/ *nds 


in the notation of Ward [10] this becomes 


_ _ 1 f W + (Rq) ■ nds j}2 f 0 + (Rp) (Rq - R) • n ds 

^ (R)- 2t ir J Rg ~2 t j J r b 3 

S w ^w 

■ by defining the conormal vector, n c =^- Bn x , n y > n z ), we also obtain 


Comparing this equation with equation ( 3 . 31 ) we see that this formula holds whatever the 
inclination of the surface S w . For superinclined surfaces, we showed in section 3.6 that 
0 and W • n take on the boundary values 0 + and W + . rt in the integrals on the right hand 
side of the equation. For subinclined panel surfaces, 0 and W • A are not independent. We 
now consider the volume V 2 which is bounded by the surfaces consisting of the portion 
of the Mach cone V on the side of S w opposite the point x, y, z and the surface T 2 of the 
lower characteristic surface from the leading edge. Since for the finite part of the integral, 
the integrals over T 2 and T vanish where the integrands are infinite, we obtain 


(0 W 2 -W 0 2 ) • n ds 


1 f W-(R 0 ).n B 2 


n J_ f * 

U 27 T J Rr 


b2 f 0 (Rp) (Rq - R) • n ds 
17 / Rb 3 
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This indicates the dependence of 0 and W • ft boundary conditions on S w for subinclined 
surfaces. This is the same surface as in equation (3.39) but ft is directed outward on the 
opposite side of S w . Adding equations (3.39) and (3.40) then yields 


^ 1 /[W+(R 0 )-W-(R 0 )] • n ds B 2 /* [> + ( R 0 ) - »-<R 0 )] (Rcr R, 

*< R) — / - “Ri -2./ r b 3 

O w b w 


Although equation (3.41 ) was derived with the point R in the upper region of 
figure 8, it is valid for a point anywhere in the flow field. We can deduce from equation 
(3.41 ) that a discontinuity in the perturbation mass flux across a surface with the per- 
turbation potential itself continuous, produces a distribution of sources of strength 

a = W+ . n - W- . n 


with the flow field described by the potential 



Similarly, a jump in the velocity potential across a surface with continuous normal mass 
flux produces a doublet sheet of strength 


M = </>+( R 0 ) -r(R 0 ) 


whose flow field is described by the potential 



f* m(Rq) (Ro-R).ftds 

< Rb 3 



Conversely, a source distribution a defined over a surface produces a jump in the normal 
component of perturbation mass flux of magnitude a and continuous potential across 
the surface while a doublet distribution produces a discontinuity in potential across the 
surface with the normal mass flux remaining continuous. 


. ftds 


(3.41) 


(3.42) 


(3.43) 
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3.8 UNIQUENESS OF SOLUTION ON SUBINCLINED SURFACES 


Proving uniqueness for certain boundary conditions applied to subinclined surfaces 
is not as straightforward as for the superinclined surfaces. Consider the following integral 
derived in Appendix 1 of Ward [10] . 


/ (u. VW. n-y V . WU. n)ds= f U . [Vv 

s v 


W- WXvxvjdv 


The volume integral becomes zero for irrotational flow in which 

v . W= vXV = 0 


We choose the surface S to consist of the three surfaces in figure 9 or 

1 . Both sides of the wing surface S w . 

2. A plane downstream of wing leading edges S3 and normal to the free stream U. 

3. The characteristic surface F\ from the leading edge of the surface S w . 

Assume U to be in the x direction. Then the integral over S3 becomes 

J* (b^u^ + + w2) dydz 

S 3 


The integral over Tj vanishes since_thejntegrand was proved by Ward to be continuous 
across characteristic surfaces and W = V = 0 ahead of the surface. Equation (3.44) 
becomes 

J' (B 2 u 2 + v2 + w 2)dydz+ J' • V+ W+ . n -y V + . W+ U • n^ds 

S3 s w 

+ / u. V-W-. n-4-V - . W". n ds = 0 
S 2 


To express the integrand in a form more suitable for considering normal mass flux and 
tangential velocity boundary conditions, we divide V and W into normal and tangential 
components. 


(3.44) 


(3.45) 
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V = (V . n) n + n X (\+. n) = (V • n) n + V t 
W=(W.n)n + nX(WXn) = (W.n)n + W t 


Then the integrand takes the form 

U . W . n - V . WU . n/2 = [u . nV . n/2 + U . V t ] W . n - W t . V t U . n/2 
and equation (3.45) becomes 

(Uq/ 2) l f (B 2 u 2 + v 2 + w 2 ) dydz + J' j [-^ U . n V + . n + U . V t + J W + . n 

S 3 s w 

-^v t +.w t +u.ft}ds {[ju.Hv-.S + u.v t -]w-.ft 

s w 



w t “ U . n 


} 


ds 


(3.46) 


The approach to proving uniqueness using the preceding integrals is similar to that used 
for the superinclined panels. The vectors V and W are chosen as the difference between two 
solutions which apply the same boundary conditions to the surface S w . If the surface 
integrals over S w are all zero, then the integral over S3 is zero. Since the integral over S3 
is the sum of squares, the integrand and the individual velocity components must also be 
zero. Since S3 is a general surface and may be moved, then u = v = w = 0 and the two 
solutions are identical. 

From the form of the integrals above, it appears that if both W . n and V t boundary 
conditions are applied to each surface, the solution is unique. However, for three dimen- 
sional configurations, we usually are interested in solutions for the flow over the exterior 
(or + side) of the surface. We then eliminate the interior flow by setting 

W“ • n = V - x n = 0 (3.47) 


Since we demonstrated in section 3.7 that the source sheet produces a jump only in the 
normal mass flux and the doublet sheet, a jump only in the tangential velocity, we have 
on the upper surface 
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W+ . n = a 
V + X n = VfiXn 


where o and fi are the source and doublet distributions on the surface, respectively. Subject 
to equation (3.47), two solutions having the same source distribution and doublet distribu- 
tion can be seen from equations (3.46) to be identical. We cannot choose both o and /i 
independently, however, since we must also have, from equation (3.47), 

n X V“ = n X v</)” = 0. 


Since this is the tangential derivative on the surface, it is identically satisfied when 

0 

or from equation (3.41), (3.42), and (3.43), 

g2 f * /z (Rq) (Rp-R).nds = j_ f nds 

"2 ST J r b 3 2tt J R b 

S w s w 


Hence, we are free to choose only one of the functions o and (jl . The most generally applied 
boundary condition is the condition that no mass flux penetrate the surface or 

W+ . n = - U • n 

which fixes the value of the source distribution, namely, 

7T A 
o = - U • n 

Uniqueness of solutions with 0 = 0 and o = -U • n boundary conditions depends only 
upon uniqueness for solutions of the integral equation (3.48) for ji. Discretization of 
equation (3.48) by a panel method leads to a set of simultaneous equations of the form 

ayjUj = b{ ij=l,2,3, ..,n(j summed) 

where i, j denotes the points at which b[ is evaluated and fx j are values of the doublet strength 
to be determined. Because of the nature of the kernel function 

B2(R n -R) . n 

R B 3 
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the ay are diagonally dominant. When the matrix is of rank n then the solution of the 
discrete values of fi are unique. When the discretization is properly carried out, then we 
have shown that the solution with the potential boundary conditions 

a = - U • n 0~ = 0 

is unique and satisfies the condition that the normal component of linear mass flux 
vanish on the exterior (+) surface. 

We shall now consider uniqueness of linearized supersonic flow solutions when 
linearized boundary conditions are applied on a plane z = 0. The quantity U • n is then 
identically zero on the surface and the surface integral in equation (3.45) becomes 


(U 0 /2) / (b2u 2 + v2 + w 2 ) dydz 

^3 


./ 


U . V + W + . n ds 


U . V W . n ds = 0 


"w 


In the same manner we let V and W represent the difference between two solutions with 
the same boundary conditions, then we see that the solution is unique if either 


U. V 


or 


W. n 


is prescribed over the surface S w . The quantity U • V is the linearized pressure and, hence, 
is the design type boundary condition while W . n = 00/3 z is the down wash boundary 
condition relating to the wing slope. From the form of the integral over S w , solutions in 
which part of the surface has design type boundary conditions and part downwash boundary 
conditions are also unique. 
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3.9 THE VELOCITY COMPONENTS NEGLECTING VORTICES AT PANEL EDGES 


The velocity components can be computed by differentiating equation (3.39). This 
approach will lead to integrals containing the doublet strength which will introduce 
vortices at panel edges when doublet strength is not continuous. Another form of the 
velocity vector which contains no integrals of the doublet strength to produce these 
vortices is obtained by using the identity established by Ward in appendix 1 of reference 
[10], i.e., 


f ( VjW 2 • n + V 2 W! . n- V! . w 2 n)ds = 


j [Vj v»W 2 + V 2 v»W] - WjX (vXV 2 )-W 2 X (vXVj)Jdv (3.49) 

*V 

If the volume is source free and the flow is irrotational, the volume integral vanishes. Let 
W 2 and V 2 be given by the source in equation (3.28) and the surface S contain a surface 
S w on which boundary conditions are to be applied, a characteristic surface from the 
leading edge of the wing surface S w , the cone from the point x, y, z and a circular cross 
section Sj cutting the Mach cone close to the singularity at xq = x, yg = y, and zq = z, as 
shown in figure 10. Since V 2 W 1 = VjW2, the last two terms of the integrand of equation 
(3.49) become 

V 2 Wj . n- Vj . w 2 n = V 2 Wj . n - "v 2 . Wjn = X (v 2 Xn). 


Now, as the distance from the singularity of at x,y,z goes to zero, we have 
f [ViW 2 . n + Wj X (v 2 Xn)] ds^ Vj J W 2 -nds + W 1 X f (v 2 Xn)ds 




J 1 


>1 


The first integral has been evaluated for the integral equation of (p and Ward [10] has shown 
that the second integral is zero, hence, 

J |V 1 W 2 . n + Wj X (v 2 Xn)]ds = 2 7 rV 1 (R) 


After substituting V] W 2 • n - Vj • W 2 n = X Vj') X W 2 , the relation for the velocity 
vector becomes, after dropping the subscript 1, 


V(R) 


R) - - 2 ~ J' [V 2 . W) + (ft X V) XW 2 ] ds 


J w 
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Substituting for V 2 and W 2 

7 <k>--£/ s - w 

s w 

If the panel is superinclined, n X V(Rq) and W (Rq ) . n take on their appropriate 
values on the surface S w . For the subinclined surface we incorporate the region on the 
opposite side of the surface lying within the Mach cone from the point x, y, z as we did 
for the velocity potential. This makes no contribution to V(R) and we obtain 


from equation (3.28) yields 


( R o) v (^) ds -#/ 

\ / s 

°w 


pxV(R 0 )]x (Rq-R) ds 


R B 3 


(3.50) 


V(R) = - 





(3.51) 


where o - W + . n - n is the source distribution on the surface and (jl = - 0” is the 

doublet distribution. The quantity n X v/i is the vorticity vector lying on the surface S w . 
Note that the doublet integral contains only derivatives of fi. This is the form for the 
velocity applied to panels when line integrals of the doublet distribution are to be 
excluded. The integral may also be written 


B_ 

2t r 



). (r 0 - rOds 
r b 3 



vt£ Rq 


-r) ♦ n 


r B 3 
(go - jO • 


ds 


vjunds 


r b 3 


(3.52) 


3.10 BOUNDARY CONDITIONS ON WINGS AND BODIES 

We have shown that unique solutions are obtained by prescribing the potential and 
the linear mass flux boundary conditions on the downstream side of surfaces inclined at 
angles to the free stream greater than the Mach angle, i.e., superinclined surfaces. When 
the flow is subsonic, the solution can be shown to be unique if either the potential or the 
normal linear mass flux is prescribed on any surface. There is no similar proof for 
subinclined surfaces in supersonic flow. When linearized boundary conditions are prescribed 
on a plane then the solution is unique when either the linearized mass flux or the linear 
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pressure are prescribed. Furthermore, the solution is also unique if mass flux boundary 
conditions are prescribed over part of the surface and linear pressure boundary conditions 
are prescribed over the remainder. 

Using either a source or a doublet distribution alone on a closed surface such as a wing 
or body, yields an interior flow with disturbances propagated along Mach lines and reflected 
from the surface and causing unrealistic variations in the exterior pressure distribution. 

By utilizing a combined source and doublet distribution on each panel, it is possible to 
eliminate the interior flow for a wing or body and to provide a smoother pressure distribu- 
tion on the exterior of the wing or body. We have shown in section 3.8 that choosing the 
source strength o according to 


xT A 

o = - U • n 


and setting 


0 = 0 


on the interior of the surface, produces the appropriate boundary conditions of no mass 
flux through the surface, 

w. n = -u. n. 


To preserve continuity of doublet strength, vortex sheets are shed from the trailing 
edges of wings. When the trailing edge is supersonic, the flow at any point on the edge is 
not influenced by any other point on the edge. The flow then changes abruptly across the 
trailing edge as in two-dimensional flow. On subsonic edges, the flow at any point on the 
edge is influenced by all points along the edge contained in the Mach cone upstream of 
that point. The flow must leave the wing smoothly and a Kutta condition is applied to all 
such edges. Special planar wake panels preserve continuity of the doublet strength with 
the wing trailing edges. The doublet strength on the panels depends only upon the 
variable normal to the flow direction. Thus, the linearized pressure jump 

AC p = - 2dp/9x = 0 


is a necessary requirement since the doublet sheet representing a shed wake must not support 
a pressure differential. 
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4.0 GENERAL DESCRIPTION OF THE SUPERSONIC PANEL METHOD 


4.1 EARLIER SIX PARAMETER SPLINE METHOD 

It was shown in section 3 that solutions to the linearized supersonic flow over a con- 
figuration can be expressed in terms of distributions of source and doublet strength on the 
configuration surface. Except for the simplest of geometries, complete analytic determina- 
tion of the source and doublet strength from the integral equation defining the boundary 
conditions is impractical. To overcome this difficulty, the panel methods were devised. 

The panel method is basically a collocation method for satisfying the boundary 
conditions at a discrete set of points in order to evaluate a corresponding set of source and 
doublet strengths. The source and doublet distributions on the panel are defined in terms 
of the unknown values of the singularities at the centers of the panel and of neighboring 
panels by a system of spline type polynomials. The resulting integrals over the panels 
representing the velocity components and potential at the control points on the panels 
can be integrated in closed form, and imposing boundary conditions yields a set of 
simultaneous equations to be solved for the singularity parameters. When these quantities 
are known, the flow field velocity and pressure can readily be computed. 

An example of a paneled wing and body combination is shown in figure 1 1. 

A sufficiently fine set of grid points is defined over the surface. By joining the grid points 
with straight lines, we form the basic set of quadrilateral panels, which in some special 
cases reduce to triangles. Generally, these quadrilaterals are not planar. One can easily 
show that the midpoints of the four sides of the quadrilateral lie on a plane which we shall 
call the average plane of the panel. In figure 12 the vector representing the side P 5 - P 5 
of the inner quadrilateral is given by 

W( ? l + P 2 )/2-(Pi +P 4 )/2 = (P 2 -P 4 )/2 


Similarly, the vector of the opposite side P 7 - Pg is 

P 7 - P 8 = (F 2 + P 3 )/2 -(P 4 + P 3 ) 12 = (P 2 - P 4 7/2 = P 6 - P 5 


We see that the sides are parallel and P 5 P 6 P 7 P 8 forms a parallelogram. Following the sub- 
sonic method of Johnson and Rubbert [3] , the earlier supersonic program utilized 
quadratic doublet and linear source distributions prescribed over the projection of the 
panel onto the average plane. Thus, the panels inducing the flow did not form a continuous 
closed surface for non-planar configurations, but contained gaps between the panel edges. 

Essentially, two boundary conditions are required at the center point of each panel 
and associated with each point are the values of the source and doublet strength to be 
determined. The six coefficients of the quadratic for the doublet distributions and the 
three coefficients for the linear source distribution on each panel, are determined by least 
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square fit with the unknown values of the doublet and source strengths at the centers of 
the panel and of neighboring panels. A system of splines then defines the doublet and 
source strength distributions over the surface of the configuration, although, they do not 
have the continuity usually associated with conventional splines. 

Because the spline fit requires surfaces with a certain degree of smoothness, complex 
configurations such as wings and bodies are divided up into networks. At junctions of 
networks, discontinuities in curvature and slope are allowed. To match doublet strength 
across network edges, additional control points at panel edges are provided. Mathematically, 
a network is a rectangular array of panels over which a set of control points, with appropriate 
boundary conditions, and source and doublet parameters are defined to provide a determinate 
system. Figures 13 and 14 schematically show control point and singularity locations for 
the source and for the doublet networks with analysis boundary conditions, respectively. 
Figure 1 5 shows the location of the control points and the singularity values for the doublet 
design network used in the single planar example described in section 12. When the flow is 
required over a known body shape, source and doublet strengths defined at these same points 
determine the source and doublet distribution over the entire network by the spline system. 

The three parameters to define the source strength and the six parameters to define 
the doublet distributions are not sufficient to ensure complete continuity of the source and 
doublet distribution at panel edges. As shown in section 6, discontinuity in the doublet 
strength induces higher order singularities in the flow than a similar discontinuity in source 
strength. The six parameter doublet distribution can be made to agree at panel corner points 
only or can be made continuous only in a least square sense across panel edges. Some 
excellent results using this spline were obtained and are discussed in section 12. However, for 
some paneling of certain configurations, the effects of the panel edge gaps and slight dis- 
continuity of doublet strength were found to build up and lead to instability in the solution. 
To distinguish this- earlier spline system from the spline system to be described in the 
following, we shall call it the six parameter spline. 

Ideally, the doublet strength at network edges should be as strictly continuous as it 
is at panel and subpanel edges. However, the doublet strength can be made almost continu- 
ous for the six parameter spline by moving the control points an infinitesimal distance away 
from the boundaries into the panel and using the infinite singularities of the aerodynamic 
influence coefficients as shown in section 6. This approach fails for supersonic leading 
edges where only finite discontinuities occur. To match the doublet strength at supersonic 
network edges, incompressible aerodynamic influence coefficients are used for all network 
edges. The inconsistency is only apparent because the role of the aerodynamic influence 
coefficient is to make the doublet strength and its conormal derivative nearly continuous. 

For greater details on the six parameter spline, the reader is referred to Johnson and 
Rubbert [3] . 
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4.2 IMPROVED NINE PARAMETER DOUBLET SPLINE AND IMPROVED PANELING 
TO OBTAIN CONTINUOUS DOUBLET STRENGTH AND GEOMETRY 


To establish a panel system which is a completely closed polyhedral surface, we divide 
each quadrilateral panel into 8 triangles. We have seen that the 4 midpoints of the sides of 
the panel lie in a plane, and the lines joining these midpoints actually form a parallelogram 
(see fig. 12). Each corner point, together with the midpoints of the lines meeting at that 
comer, form a triangle which has a common edge with the parallelogram. Thus, the con- 
ventional quadrilateral panel can be divided into 5 planar surfaces whose edges are contiguous 
and which are contiguous with the edges of adjoining panels. When the center parallelogram 
is divided into 4 triangles, the doublet can be made continuous across all panels and sub- 
panel edges by using the nine points in figure 12 to define a separate quadratic distribution 
in each subpanel. Because of this feature, we shall call this splinal system the nine parameter 
spline. Since the center of the panel where values of the source and doublet strength are 
defined is the common vertex of the four panels, the control point at which the boundary 
conditions are applied must be displaced from the center to avoid numerical difficulties 
with the panel influence coefficients which are indeterminate at panel comer points where 
s = s m = z : =R = 0 (see, for example, equations (A1 5) and (A1 6)). Since the nine parameter 
spline renders the doublet strength continuous and there are no gaps in the edges, the 
aerodynamic influence coefficients can be simplified by eliminating all line integrals of the 
doublet strength around panel edges as these integrals will cancel in the summing of the 
aerodynamic influence coefficients from the adjacent panels. For greater details in the 
construction of the spline, refer to appendix C. 

The nine parameter spline insures continuity of doublet strength everywhere within a 
network. In the latest version of the program, the doublet is made strictly continuous also 
across network edges, by appropriate boundary conditions at the control points on the 
network edges (see fig. 14). Consider two networks, the second lying downstream of the 
first. On the downstream edge of the first network, regular downwash boundary condi- 
tons are applied to determine the values of the doublet strength. These values of the 
doublet strength are then used to define the doublet strength at these control points common 
with the adjoining downstream network. By means of the spline system, the doublet 
strengths are defined at the panel corners along the network edge in terms of the values at 
the network edge control points, and are restricted to have the same values for both net- 
works. Since there are three points on each panel at which the doublet strength is defined 
along the network edge and a cubic spline is used, the doublet strength is then made con- 
tinuous on all points of the common boundary between the two networks. 

Since some of the example flows discussed in section 1 2 were calculated using less 
strict matching of the doublet strength, it is worthwhile, therefore, to discuss this indirect 
method of nearly matching the doublet strength across network boundaries. In section 6, 
an explanation of how the singularity in velocity from the line integral of doublet strength 
across network edges is presented. Since these line integrals have been eliminated from the 
aerodynamic influence coefficients used in the 9 parameter spline, a simulated incompressible 
vortex was added to the conventional mass flux boundary conditions. To achieve a match 
of the doublet strength across network edges, the boundary points were moved an infinitesimal 
distance into the network and away from the panel boundary where the influence coefficients 
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become infinite like 7 /r where r is the distance from the edge and 7 = Am is the vorticity 
magnitude. Thus, at two points on opposite sides of the network boundary we have 


(w 1 +U).ni+(Ml-M 2 )/ri =0 

(W 2 + D) . n2 + (M2 -J“l)/r2 = 0 (4.1) 


where rj and ^ are distances from the points 1 and 2 , respectively, to the network edges. 
These equations may be combined to give 

M] - ju 2 = (Wj + U) • njr] - (W 2 + U) • n2r2 ~ 0 


and 


r l (^1 + U)*n + r2 (W 2 + U) • n2 = 0 


(4.2) 


Hence, the doublet strength is nearly matched and mass flux boundary conditions are 
satisfied in a weighted average. 

4.3 COMBINED SOURCE AND DOUBLET PANELS 

On aircraft bodies and wings without lift, it may be possible to find the linearized 
solutions with sources alone on the surface. For such closed configurations there results 
both exterior and interior flows. When the flow is supersonic the interior flow may cause 
build-up of disturbances which propagate along Mach lines and are reflected from the 
surface in the downstream direction as shown in figure 16. These waves may produce large 
pressure disturbances on the exterior of the surface. Such interference waves can be 
eliminated by combined source and doublet paneling on the configuration. 


In section 3.8 we showed that when the source strength at the control points is set 
equal to -U • ft, i.e. 


a = - U • 


A 

n 


(4.3) 


and the boundary conditions 


0 = 0 


(4.4) 


are prescribed at points on the inner boundary, the boundary conditions 

W . n = - U • n (4.5) 

result on the outside surface. Furthermore, the interior perturbation flow is zero, 
therefore eliminating the interior waves and the resulting disturbances on the exterior 
surface. 
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The application of boundary conditions in equations (4.3) and (4.4) on the surface of 
a flow-through nacelle will only work on an isolated nacelle in supersonic flow where there 
is no incident perturbation flow, for which 0 is automatically zero on the inlet surface. 
When a nacelle is situated in the region of influence of a wing, the inlet can be closed by 
means of superinclined panels, which absorb or cancel the incident flow. Boundary 
conditions can only be applied to the downstream side of the superinclined surface; and to 
eliminate the incident perturbation flow, these boundary conditions, as explained in 
section 3.6, are 


0- w.ft =0 


It is also possible to represent a jet exhaust by terminating a nacelle with a super- 
inclined network. For this case the normal mass flux would be prescribed, namely, 

w.ft=f(€,a) 

The jet exhaust cannot be represented exactly, since the issuing jet has the same total 
head as the free stream flow because it satisfies the same differential equation. A convenient 
second boundary condition on the exit surface is 0 = 0; and to match doublet strength a 
trailing doublet wake must be included. The shape of the sheet is not known but it may 
be approximated by a cylinder. 

Thin wings can also be represented by sources and doublets distributed over the planar 
planform. The thickness distribution is represented by sources, while camber shape (or 
angle of attack) is defined by doublets. Examples of cambered wings and flat plates at 
angles of attack and of symmetrical thick wings are presented in section 12 along with 
comparison with exact linearized solutions. 

4.4 EFFICIENT COMPUTATION OF THE INFLUENCE COEFFICIENTS 

The influence coefficients were derived in such a way that they give valid results when 
all or part of the panel edge lies outside the zone of influence of the control point. However, 
since their computation is fairly costly, considerable economy in computation is achieved 
by using a subroutine to test whether a panel is in the zone of the given control point and 
eliminating computation of the influence coefficients for those panels lying outside. The 
first test is to find out if any comer of the panel is within the upstream Mach cone from the 
control point by computing the square of the hyperbolic distance from each comer. If it 
is positive for any corner, then computation of the influence coefficients is carried out. 

If it is negative for all corners, then a further test is made to see whether any edges cuts the 
Mach cone. If so, then the computation of the panel influence coefficients is carried out. 

If not, then the program proceeds to the next panel without entering the influence coefficient 
subroutine, unless the panel happens to be in a superinclined network. For superinclined 
panels, a further test is required to find whether the intersection of the Mach cone with the 
panel lies entirely inside the panel. 
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A further economy in computation is achieved by recognizing that some terms of the 
influence coefficients are linear combinations of a smaller set. Also, since the terms such 
as equations (A15) and (A 16) are evaluated at the end points of the edges, the inverse 
tangents and logarithms from the two endpoints of each edge can be combined, requiring 
one half the necessary entries into the arc tangent and logarithm subroutines, which are 
costly in computing time. 

For panels far away from the control point and lying well within the Mach cone, far 
field expansions of the aerodynamic influence coefficients may be used. Since these are 
algebraic in form they require less computing time than the near field aerodynamic influence 
coefficients. 
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5.0 AERODYNAMIC INFLUENCE COEFFICIENTS FOR 
SUBINCLINED DOUBLET PANELS 


To represent the aircraft configuration with sufficient fidelity by means of the simple 
panels described in the foregoing section, it must be divided into panels sufficiently small 
that the curvature of the panel is small in relation to the panel’s smallest dimension. The 
panel shape can then be represented by a simple second degree polynomial. Linear distri- 
butions of source strength and quadratic distribution of doublet strength are utilized 
enabling the contribution to the velocity potential from each panel to be integrated in 
closed form. 

To obtain the simplest form for the influence coefficients, we scale the variables 
according to 


(x,y,z) = (x c ,By c ,Bz c ) 

where x c , y c , z c designates the compressible coordinate system, so named because x c is in 
the direction of the free stream velocity. This transformation eliminates the factor B from 
the hyperbolic distance Rg. For the six parameter spline, we represent the subinclined 
panel as a small quadratic departure from the plane z = 0 and integration is performed over 
the projection on the plane z = 0. For the nine parameter spline, the subpanels themselves 
are planar and are assumed to lie in the plane z = 0. Choosing the subinclined panel as the 
z = 0 plane also simplifies the conormal derivative to 3/dz, further simplifying the integral 
for the doublet. For subpanels or panels not lying in the z = 0 plane, the influence 
coefficients are transformed by coordinate transformations which are discussed in section 10. 

The integration over the panel may be accomplished by combining the integrations of 
regions bounded by a panel edge, two lines of constant y and the intersection of the Mach 
cone on the plane z = 0 as shown in Figure 17. This can be seen by considering figure 1 8 
for a panel which lies wholly within the Mach cone. The integration over the area between 
lines AB and BC are added together. Since CD and DA are taken in the reverse order, the 
integration over the area between these two lines are subtracted from the areas under AB 
and BC. The area remainding is the desired area of the quadrilateral ABCD over which 
the integration is required. The integration over the region in figure 17, yields a quantity 
to be evaluated at the endpoints of the panel edge. The induced flow from the panel 
is then computed by evaluating the function at the endpoints of each side proceeding 
around the panel in a counter-clockwise direction. 

Let the panel surface be defined by 

f ( x o>yO’ z o) = z o - f ( x o>yo) = o (5.i) 

Then to the first order in f , the normal to the surface becomes 

n ZQ = 1 + o(f^), n X Q = - f X Q + 0(^2^ + o(f2) (5.2) 

where r XQ = 9r/9x 0 and f yQ = 9?/9 y() . 
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We shall consider the doublet term first. Thus, with i// = 1/R S where 
R s = \/( x_x o) 2 " (y-yo) 2 -( z - z o) 2 > e q uation (3-43) becomes 

i**l**«“ &/'/(% + fx 0 ^-?yo3^)"< xo ' yo )‘ iX| > dyo 
s w s w 

where ds = dxQdyo/cosft = dxQdyg + 0(f 2 ) , 3/3n c = (- n x , n y , n z ) • v , n = (n x , n y , n z ) 
is the normal of the wing or body surface, and 12 is the angle of projection of the element 
of the area on the average plane. The sign of equation (3.43) has been changed to conform 
to the general practice of defining the body and wing normals outward into the fluid, while 
the development of the general theory in section 3 is based on an outward normal from the 
fluid region. 

Since the finite part contains no terms from differentiation of the limits of integration, 
we may move the differentiation outside of the integral sign by noting that 


3 _ 3 \p 3 \p _ 3 djjj _ _ d\p 

3xq dP 3yg 3y 3 3zq 3z 


Therefore, for the velocity potential we obtain 


^ 2tt 3z 




m(xq, yo)dx Q dyQ 
R s 


JL JL 

27 r 3x 




m(x q, yo) ^xn dx Q d yo 

R s 


with 


+ 


j_ JL 

27t 3y 



siyo) 


C ( y o) n ( XQ , y Q ) f y dx 0 dy 0 


Rc 


(5.4) 


(5.5) 


R S =y( x - x o) 2 - (y-vo) 2 - ( z ~0 2 - 


The finite part symbol is removed when integration is performed before differentiation 
since the integral is then made convergent. Here 

? = ? ( x o> yo) and x o = f (yo) 
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is the equation of the projection on the plane z = 0 of the intersection of the Mach cone 
with the panel surface. The region of integration is shown in figure 17. We simplify the 
Mach cone singularity by introducing a new variable for xq. Let 


x o = £(vo) - 1 = £o - 1 


and assume a straight edge 


x o = x i + (yo-yi) / m 

where m = ^y 2 - yj) /'(x 2 - xj) is the slope of the panel edge. By this transformation 
t = 0 represents the Mach cone boundary. The hyperbolic distance becomes 


R S = y( x -£o + 0 2 - (y-yo'' 2 - [ z -f ($o - yo)] 2 

Expanding in the form 

f (£o “ L yo) _ ? (?o> yo) ~ yo) + ^fxx/^ 

and using the condition that 

/Ho) 2 - (y-yo) 2 - [ z -f (So. yo)] 2 = 0 


yields for the hyperbohc distance 


yt[aj t+b, +o(r 2 X 


where 

a l = 1 + ( z ”fo)^xx ” 1 + z fxx + 

b l = 2 [( x ~£o) - z ?0x] 
fo = f(£o>yo) 


(5.6) 


(5.7) 


(5.8) 
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To find the equation for the Mach line boundary we consider 


( x -£o) = v/(y-yo) 2 + ( z- ?o) 2 = r - z ?o/ r + °fr 2 ) 

or £o = x-r + z f()/ r + 0 (f 2 ) 


(5.9) 


where r= .^/(y-yo) 2 + z2 - Th en ^ (?0= yo) ( x-r > y o) + ®(f 2 )- 
Finally retaining only first order terms in f yields for the potential 


<t> = 


2tt 



(yo >) /m+e i ^ ^x-r-t+e], y 0 )dtdy 0 

(a 1 t+b ! ) 



r ^ y ° yi )' /m r x ( x ~r~t, ypV(x-r-t, yp^dtdyp 

x /t(t+2r) 


where 




x_x i - r ~(yo-y i) / m 


?y (x-r-t. yp M(x-r-t, y 0 )dtdy 0 
V / t(t+2r) 


(5.10) 


ej=zf/r e 2 = z $xx e 3 =z Cx+f/ r ) 

a l = l+e 2 b] = 2(r-e 3 ) (5.11) 


The integration of equation 5.10 is carried out in appendix B retaining only the first 
order terms in the epsilons and f . It was found that the contributions to the velocity by 
the first order terms have higher order infinite singularities near the Mach cones emanating 
from panel corners and hence, are not a true correction for panel curvature. Consequently, 
the influence coefficients for the flat panels only are used. Setting f = 0 in equation 5.10 
yields 


, = ii 

^ 27 r dz 



-r-(yo-yi) /rn 


ju(x-r-t, y 0 )dtdy Q 


(5.12) 
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The integral may be simplified by introducing the following variables: 


x m = x - X 1 “ (y - yi)/ m > s = y-y 0 
s m = x m + s/m 


( 5 . 13 ) 


Scaling the t variable by r eliminates differentiation of the singular part of the integrand 
allowing differentiation to be performed before integration. We obtain for the perturbation 
potential 


* = - 


JL 

2tt 



p(x-r-rt, y-s)dtds 


( 5 . 14 ) 


Differentiation leads to 



y-yi 


s m M(x-s m , y-s) ds 

r2 \/ s m 2 - 1-2 



ju x (x-r-t 5 y-s) (r+t)dtds 
r 2 v/t(2r+t) 


The integration with respect to t may be performed by expanding p x in powers of t using 



t n dt 

v 't(2r+t) 


Hence 


= -— f 

2tt J 


y V1 ( s m h(x-s m ,y-s) PxQq y-s) (s m 2 -r 2 ) 


y— y 2 


r 2 R 


r 2 R 


. , s m ( s m 2 -r 2 ) 

- M XX IQ + r 2 R 


/2 > ds 


( 5 . 15 ) 
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where we have eliminated I 2 and 1 2 by the recursion formulas in equation (A8). Finally, 0 
is obtained by expanding ju about x, y and performing the integration with respect to s. 
Using the functions defined by equations (A4), (A5), (A10), and (A12) leads to 


<t> = - 


z 

27 r 



-Mx w 0-My Qml + M x x( w l/ m + z2 Qm0)/ 2 


+ M X y w l 


+ Myy ( w m0 z2 Q m o)/2J 


(5.16) 


The variables zQ m Q and Q m i may be simplified by discarding those terms which depend 
only upon the corner point and therefore cancel when integrating around a panel corner. 

The simpler relation takes the form 

0 = - 2 “ [juQi -M x zwq - Myzwo/m + m X x z (w]/m+zQi)/2 

+ M X y zw l +Myy z(w m o-zQi)/2] (5.17) 


where zQ m o has been replaced by the function Qj defined in appendix A by equation (A 19) 
and Q m j by w Q /m in equation (A18). We have also discarded the term sIq from the func- 
tion S 0 since it cancels when integration is performed around a panel comer. 
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6.0 PANEL EDGE SINGULARITIES FOR DOUBLET PANELS 


To investigate the singularities which result from discontinuities in doublet strength 
and its derivatives along panel edges, we consider the potential from the doublet distribu- 
tion in equation (5.17). The perturbation velocity vector then becomes 


v 0 = - 27 [m vQl + QiVju - (*i x + My/m) v (zwq)- zwq v (ji x + n y /m ) 

+ (n xx /2m + n xy + n yy /2m) v(zw])+ M yy v (zx m WQ - z 2 Qj) / 2I 


( 6 . 1 ) 


The contribution to the velocity from a specific panel at a given point (x, y, z) on any 
panel is found by evaluating the functions Qj, wq> and wj and their gradients at the endpoints 
of each panel edge going around the panel in a counterclockwise direction. For interior 
points in a network, the induced velocity from a single edge has contributions from the two 
adjacent panels. If the representations of n on the two panels are equal at the panel edge 
then the term involving m in equation (6.1) cancels when the reference point approaches 
the common boundary. However, a jump in the doublet distribution or in its gradient at a 
panel edge introduces singularities in the velocity at the panel edge. To investigate these 
singularities, we study the behavior of vQj, v(zwq) , and v (zwj) as the point (x, y, 0) on 
the panel approaches a point on the panel edge. From equation (A33) for z = 0, we obtain 


vQi = 



where s m , s correspond to x-xj, y-yj or x-X2> y-Y2 f° r ea °h s ^ e * 

Now x m may be interpreted as the coordinate measuring distance from the panel edge 
which is defined by x m = 0. The oblique coordinate is proportional to y m = s - s m /m. Thus 

vQ! = (0, 0, y m /x m v/sm 2 - s 2 ) (6.2) 


On z ~ 0, 


V(ZW() = (0, 0, Wi) 


(6.3) 


For points close to subsonic panel edges, we have shown in equation (A3 1) that 


w 0 ~ ~ / == ■ !°g x m 

VI -m 2 


(6.4) 
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and from equation (A9) 


m^ x m log x m + m^R 

7(1 -m2) 3 1-m 2 


(6.5) 


Note that only 80/dz is singular. Hence, substituting equations (6.2) through (6.5) into 
equation (6.1 ) yields for points near the subsonic panel edge on the panel 


0 Z 


27T 


{ M(x,y )y m /x m Vs m 2 - s 2 - (M x m + M y ) lo S x m/ \A ~ m2 + 0(1)} (6 - 6) 


The oblique coordinates associated with the panel edge represented by x m = s m - s/m = 0 
are defined by 


x' = (- ms m + s) / s/\ - m 2 
y' = (- ms + s m ) / 7l -m 2 


With these coordinates, equation (6.6) may be written 

</> z =- 2 ~ j^y'/x' 7 s m 2 ~ s 2 -OW9x')log x' + 0(1)} 

A discontinuity in fi is seen to introduce a singularity of 1/x' as x' goes to zero near 
the panel edge x m = x' = 0, while a discontinuity in the conormal derivative 3/i/8x' intro- 
duces a logarithmic singularity. As the panel edge approaches a Mach line in direction, the 
conormal derivative of n becomes the tangential derivative, and the 1/x' singularity reduces 

to lA/xVsince y' = x\ The continuity of ii then eliminates both singularities when the 
panel edge is a Mach line. 


For subsonic panel edges these singularities furnish a means of enforcing continuity 
of ju and of its conormal derivative at network edges by employing ordinary downwash 
boundary conditions near the edge of the network. Consider the boundary conditions at 
two very close points on opposite sides of the network edge. We assume linearized boundary 
conditions for simplicity in demonstrating the effect of the singularity. At point 1 , we have 


a</>/9z = - ^ j [mi (x'l.y'i) -n ( x 'i . y'i)] y'i/ x 'i R i 

+ [mix' (x'l.y'i) -M2x' ( x 'l .y'l)] logx'i +0(1)) = (9z/dx) x 'j (6.8) 
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or 


[mi (x'i »y'r) - m 2 ( x i > y'i )] y'i / R i 

= [mix' ( X 1 .y'i ) - M2 x'j ( x 'l > y'i )] x l 1°8 x 1 -2tt x'j (3z/3x) x ', (6.9) 

where 3z/3x is the local streamwise slope of the wing and the subscripts 1 and 2 denote {i 
as defined by the parameters associated with the two adjacent panels having the edge in 
common. The right hand side of equation (6.9) is negligible for x' very small and we then 
approximate the condition 

fi2 ( 6 . 10 ) 


at one point on the panel edge. 


At the second point on the other side of the panel edge, the continuity of (i by equation 
(6.9) is satisfied to the order of x m log x^. Then the boundary condition at point 2 is 
dominated by 

£ [*l-M2]“-2»(H) 2 /lo g x' 


from which we have 


dn\/dx ' — 3jU2/9x' 

thereby enforcing near continuity of the conormal derivative of fi. 


We now consider 30/3z near a supersonic panel edge. For this case, derivatives of 
WQ, W] , and Qj vanish ahead of the Mach cone near the panel edge (they actually become 
imaginary); but wq, wj, Qj take on constant values inside the Mach wedge shown in 
figure 19. From equations (6.1), (A27), and (A30), we see that on the panel near the 
supersonic panel edge 




2^ L“ (^x + My/m) w 0 + ( Mxx/ 2m + Mxy + M yy /2m) w 0 


Myy 


-J i 
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or 


$7 


— | M x ( x m/^) + ^Mxy + ^Myy) V (1-X?) “ Myyj | 


For the supersonic leading edge the coordinates along the panel edge and oblique or 
conormal to the panel edge are 


x ' = ( s m “ Xs ) / = x m / 

y' = (- Xs m + s) / >/l - X 2 


Hence 


. _ 1 

^"-4 


+ O(x') 


For a supersonic leading edge, the jump in the conormal derivative of the doublet 
distribution introduces a jump in the velocity conormal to the plane of the panel. 

In addition to the singularities on the panels near the panel edges, there are infinite 
square root singularities which occur on the Mach cones^emanating from panel comers out 
into the flow field whenever the doublet distribution is not continuous at panel edges. 
From section A 6 of appendix A, it is easily seen that the gradients of the functions have 
the hyperbolic distance R in the denominator and hence, become infinite on the Mach 
cone R = 0. If (jl is continuous these singularities are cancelled, providing the edges of 
two adjacent panels coincide. Although this is not readily apparent from the form of 
equations (6.1), it is easily shown in section 7 where the line integrals of the doublet 
distribution along panel edges are discarded in determining the potential and velocity 
components. 


48 



7.0 INFLUENCE COEFFICIENTS FOR SUBINCLINED DOUBLET 
PANELS NEGLECTING THE EDGE VORTEX 


It has been shown that jumps in the doublet distribution across panel edges introduce 
singularities in the velocity which propagate into the region away from the panel. The 
original six parameter spline provides a doublet distribution which is almost continuous 
across panel edges. Since the discontinuity is small, it would appear that the difficulty 
could be removed by arbitrarily eliminating the line vortex integrals from the velocity 
aerodynamic influence coefficients. This was quite successful in subsonic flow, but it was 
found that for supersonic flow even small discontinuities in the doublet at panel edges 
produced inaccuracies when the line vortex was omitted. For the nine parameter spline, 
the doublet distribution fi is continuous and we can therefore simplify the aerodynamic 
influence coefficients by eliminating all terms involving line integrals of doublet distribution 
along panel edges. In the influence coefficients previously derived, these integrals are included 
but introduce no difficulty if the doublet strength is continuous along panel edges. In the 
following we derive influence coefficients for the flat panel in which line integrals along 
panel boundaries of the doublet distribution are eliminated. Consider the doublet distribution 
from one panel and the resulting velocity potential from one edge. This is given by equation 
(5.14) or 




/i(x-r-rt,y-s)dtd s 
V / t(t+2) 


(7.1) 


Performing the differentiation yields 


r y-y\ r 

-*/ / 


s m“ r 


(r+t)/i x (x-r-t,y-s)dtds 


y_y 2 


r 2 ,/t(t+2r) 



y— y 2 


M (x-s m , y-s) ds 
r 2 %/s m 2 -r 2 


The second term is’a line integral of the doublet strength along the edge. Hence, we discard 
it and obtain 


0 = - 


2tt 


y-yi 

/ 

y-y 2 


I 


>m 


(r+t)/i x (x-r-t,y-s) dtds 
r 2 v /t(t+2r) 


(7.2) 
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The velocity potential is obtained by expanding ju x in powers of r + 1. Thus we obtain 

rY Y2 ^* Sm r [p x (x,y-s)(r+t) - M xx (r+t) 2 ] dtds 
y-y 2 *0 r 2 y/Ut+ 2) 


* = ~T* 


Z 

2tT 


/y-y2 

[ M x (x,y-s) v / s m 2_ r 2/ r 2 _ Mxx ( r 2 lQ + Sm R) / 2r 2] ds 


y-yi 

where we made use of the recursion formulas for I n in equation (A8). After integration with 
respect to s, 0 becomes 

0 = ~ 27 t [^x^O “ ^xy^l ” ( z ^0 + 

= - 2 ^[ M x p 0 -M X yPi -^ xx (zx m w 0 -z 2 Q I + x m P 0 + P 1 /m)] (7 - 3) 

where P n = z ^Rs n ds/r 2 . 

The relation for 0 in equation (7*3) ultimately was not used since it involves the P n 
functions which are not required for the velocity components. The expression for 0 in 
equation (5.17) is valid for continuous ju, introduces no infinite singularities for discontinuous 
M, and involves the same functions as the velocity components for continuous fi. Hence, it 
is more economical to use equation (5.17). 


For the perturbation velocity component / ul we have 


*x = -2^ 


•y yj /• s m _r 


h xx 


/ 7 


(r+t)dtds 


y_y 2 


/ 


y Yl s mPx (x-s m , y-s) ds 


r 2 Vt(t+2r) y ~y 2 


2 7 s m 2 “ r 2 


Integration yields 


0x = " £ { MxQ m 0 - Pxy w o/ m " ^xx^o} 
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or 

~ " 2? { M x Ql ~ Mxyzwo/m - M X x zw 0 f ( 7 - 4 ) 

For the y component of the perturbation velocity, we have 



The last term can be discarded since the coefficients of p x , p x y , and p xx from this 
term depend only upon the endpoints and not on the panel edge slope m, and therefore, 
will cancel in the integration around a panel. The quantity p x in the first integral can be 

replaced by p y and the derivative along the edge p (x-s m , y-s) which will cancel when 
p is continuous on panel edges. We thus obtain 



Integration yields 

^y = ~ 27 T | ^yQmO “ Mxy QmmO “ ^yyQm j + Mxy [QmmO “ w o] 
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or 


_ ~ 2 It { M yQl " y zw 0 ” M yy zwo/m} 


(7.6) 


To find (j> z we consider the potential in the original form and obtain 


= _L 32 /* y ' yi f 

^ 3z^ J J 

y-Y2 o 


M(x-r~t,y-s)dtd s 
y/t( t+2r) 


Since the operator 3^/3^ equals the operator 3^/3x^ - 3^/3y^, we may write 


0 Z 


.1 Ail 

2^ \3x^ 



M(x-r-t,y~s)dtds 

>/t(t+2r) 


The first integral becomes 

jl r yi r 

3x2 y y 


m 


y-y2 


o 


n(x-r-t,y-s)dtds 

>/t(t+2r) 


_3 

3x 



y-y 2 


m (x-s m , y-s)ds 
V^m 2 ~ 1-2 


* / 

y-y 2 0 


M x (x-r-t, y -s)dtds 
v/t(t+2r) 



y— y 2 


M x ( x ~ s m> y~s) ds 
x/ s m 2 “ f2 


+ i“xx s 0 ~ ^x w 0 - ^x y wj - M xx ( w m0 " S o) 


Since s 0 sI 0 + x m w 0 QmO z m0 2 and sIq may be discarded as it depends onl y upon the 
endpoints, we get for the integral 


- p x wq - M xy w j - p xx (w j/m + zQj ) 


(7.7) 
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The second term becomes 


0 r y-y\ /• s m -_r 

f I u(x-r-U v-s)dtds 

3y 2 ' J v/t(t+2r) 

y-y 2 0 v 



y-yi 

u(x-r-t, y-s)dt 

y-y 2 




M v (x-r-t, y-s)dtds 
v/t(t+2r) 


J_ /* y y * ** (y~ s m> y-s)ds ) 

m y 4 2 I 


The last integral on the right hand side of the foregoing equation vanishes because of 
continuity of ji on the panel edge and because its coefficients of ju x , etc. depend only upon 
the panel endpoints and hence, cancel in the integration around the panel comer. The first 
integral becomes 


JL 

0y 



( yo~y\ V 
M V x i — ~ -yp/dyo 


^( x_x i +x iP~) 2 ~ (y-yo)' 


This integral vanishes because of continuity of [i. The remaining integral yields 

y_y l My(x— s m , y-s)ds 

y- y 2 


i f 




MyySQ 


(7.8) 


after discarding terms depending only upon endpoints. Finally, <t> z becomes, after combining 
equations (7.7) and (7.8), 


0z = 27 l^x w 0~^xy w l “ Mxx(w]/m+zQ! ) 

+ My wo/m -M xy w m o/m-Myy(SQ + w^m)}- 


Since w m o/m = x m wo/m + wj/m 2 and from the recursion relation in equation (A9). 

x m w o/rn = R - (l - m 2 )wi/m 2 

we see that w m o/m may be replaced by wi since the term R may be discarded. Thus, we 
obtain finally 

= i l MxW o + ^y w o/ m - 2 ^xy w 1 

-M X x(w]/m + zQj) + Myy(zQi- w m0 )} (7.9) 
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8.0 AERODYNAMIC INFLUENCE COEFFICIENTS 
FOR SUBINCLINED SOURCE PANELS 


From equation (3.42) the velocity potential for the source distribution is given by 


d> = 


-iJf 


g ( x 0> yp) dxpdyo 


V ( x - x o ) 2 - (y-yo ) 2 - (z ~^ 2 


( 8 . 1 ) 


where z Q = f (xq, yo) is the panel shape, assumed to depart only slightly from the plane 
z Q = 0. The sign of equation (8.1) was chosen opposite from equation (3.42) because we shall 
use a normal outward from the wing or body into the fluid instead of outward from the fluid 
region as used in the theory of section 3. Introducing the variable t as in the doublet potential 
for integration between an edge and the Mach cone in figure 17, we obtain the integral 


0 “ 


± If 

y i o 


y 2 f x-x i -r- ( y o - y i) / m+e 1 o(x . r . t+ei>yo ) dtdyQ 


\/t( a i t + b 0 


( 8 . 2 ) 


where aj, bj, and e\ are given following equations (5.8). This is easily seen by comparison 
of the first integral of equations (5.5) and (5.10), respectively. 

The integration of the velocity potential including all first order terms in panel curvature, 
is described in appendix B. Since the terms related to the curvature introduce higher order 
singularities, we utilize only flat panels in the supersonic panel method. Setting f = 0 in 
equation (8.2), we obtain 



cr(x-r-t,y-s) dtds/ y/t(t+2v) 


(8.3) 


Since a is linear, we have 


a (x-r-t, y-s) = o (x-r, y-s) - a x t 


and 0 becomes, after integration with respect to t, 



y_y 2 


[a(x-r, y-s)Io - ff x Ii] 


ds 


where I n is defined in equations (A8) and (A 14) in appendix A. 
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Expanding a further yields 


1 /* y_yi 

= 2n J [ff(x,y-s)I 0 - ct x r] ds 


y-Y2 


(8.4) 


where we made use of the recursion formula for Ij in equation (A8) of appendix A. With 
the functions defined in equation (A4) of appendix A, the velocity potential finally takes 
the form 


0 = 2 ^ [tf(x,y)S 0 - ct x R 0 - a y Si] (8.5) 

where the functions S n and R n are defined in equations (A1 1) and (A12) m appendix A. 

The velocity potential can be simplified by dropping all terms in Sq, Sj which depend 
only upon the endpoints and, hence, cancel in integrating around a panel corner. This leads 
to 


0 = 27 { a ( x ’y) ( x m w 0 - zQl)- a x [ ( x m 2 - z2 )wo 
+ x m w i/m]/2-(7y( x r n w i -z 2 w 0 /m)/2j- 


( 8 . 6 ) 


The velocity components are found by taking the gradient of 0 in the form of equation 
(8.5), namely 


^0 = 27 I v ^Sq + ovSq - o x vR 0 - o y vSj } 


(8.7) 


From equation (A24), 


vS 0 = (w 0 , I 0 (r,s m -r) - w 0 /m, - zQ m0 ) 


This may be simplified by discarding terms depending only upon s, s m , and z.and using 
equation (A20). Hence, 


^Sq = (w 0 , - WQ/m, - Qi) 
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Similarly, with the aid of equations (A24) and (A18), we have 


vSj = (wj.-wj/m.-zwo/m) 
vR 0 = ( w m0> - w m0/ m > " zw 0) 


This may be further simplified by noting that 

w m0/ m = x m w C)/ m + wj/m 2 


and, from equation (A9), 


w mo/ m = R - (l-m2)w]/m2 


Discarding R, then yields 


vR 0= ( w m0> - w l= - zw o) 


Finally, we obtain the following perturbation velocity vector contributed by one panel 
edge when evaluated at the endpoints s = y - y \ and s = y - y2- 


7 <f> = 2~ j<J(x,y) C w 0> " w 0/m» - Ql) - (zQi + w i /m, - w i , - zwq) 
- CT y (w j , zQ r - w m0 , - zw 0 /m) | 


On the plane z = 0 which represents the panel, we have 

d(j)/dz = - a(x,y)Qj/27r = ± a(x,y)/2 
when the point x, y lies between the ends of the panel edge, and 

30/6 z = 0 

outside of this region. This is the basis of applying the analysis boundary conditions for the 
planar wings. The difference between the upper and lower slopes is seen to be 

30 + /3z - 30~/3 z = a(x,y) 

and the source panel represents the thickness distribution for the planar wing. 
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9.0 AERODYNAMIC INFLUENCE COEFFICIENTS FOR 

SUPERINCLINED PANELS 

9.1 VELOCITY POTENTIAL FOR SOURCE AND DOUBLET PANELS 

We have seen that to represent solid surfaces, doublet and source panels must be inclined 
at angles less than the Mach angle. Singularity surfaces for panels at angles with respect to 
the free stream greater than the Mach angle do exist theoretically, however. These are useful 
for closing inlets of jet engines and representing boundary conditions at the exit of a jet 
engine. The intersection of the Mach cone with superinclined planes are circles and ellipses 
and the integrals of the aerodynamic influence coefficients over a panel must be computed in 
a somewhat different manner than for the subinclined panels. We shall assume that the 
panel is aligned normal to the free stream direction. All other inclinations can be treated by 
the use of a rotation of coordinates and an oblique transformation, which is discussed in 
Section 10. 

The region of integration of the influence coefficients has several special cases. Figure 
20a shows the intersection of the Mach cone lying wholly within the panel. No edge or comer 
is involved in the computation and the velocity potential and velocity take on their simplest 
form. Figure 20b shows a Mach cone cutting one side of the panel, but all comers lie outside 
the domain of influence. Figure 20c shows one corner of a panel lying inside the intersection 
circle of the Mach cone. The potential contributed by any panel can be computed by adding 
and subtracting the integration over the three basic areas illustrated in figure 20. Consider 
figure 20d. The integration over the panel may be found by subtracting, from the basic 
integration over the circle, the integral over the area bounded by an arc and by a line cutting 
across the circle as shown by the three sides as extended in figure 20d. The triangular region 
formed by extending the two sides of the triangle must then be added since it is subtracted 
twice. 

Another approach is to use polar coordinates in the integrals. By this means it is 
possible to reduce the integration over a panel to a contour integral over the edge of the area 
common to both the Mach cone and the panel which we shall designate by E (see fig. 20). 

To evaluate the potential for the source and doublet distribution over a panel, we consider 
the fundamental source potential from which we can obtain the doublet potential and the 
velocity components by differentiation. This may be expressed in the form 

2tt0 s = f df • R s 2 = x 2 - (77-y) 2 - (f-z) 2 (9. 1 ) 


Now the singularity distribution fi expanded about the control point x, y, z, becomes 

m (v£) = M(y,z) + n y (v-y) + m z (?- z ) + Myyfa -y) 2 /2 

+ MyzOi-yXf-z) + m zz (?-z) 2 /2 
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This may be written in vector notation as 


mO?,?) = M + v/i . R] + Rj . v vju . Rj/2 (9.2) 

where V = 3/3y, 6 / 6 z operates only on iu and R\ is the vector (r? - y), (f - z). Substituting 
equation (9.2) into equation (9.1) yields 

2n<t> s = f [m + v/i . Rj + Rj .v v/jl . Rj/2] d? 7 df/R s ( 9 . 3 ) 


now R]/R s = -V'R S where V' = 6 / 677 , 3/3?- Substituting into equation (9.3) yields 
27T0 S = j [m/R s - vm • v 'R s -~Rj .v vju . v 'R s ] d 77 df 

•'x 


j* [m/R s + v - v' • (vjuR s ) - * v' . (Rj .v vmRs)]^ 7 ^? 


We now apply the divergence theorem to the last two integrals and obtain 


27T0 S - C [m/R s + V • v uR/ 2 ] drjdf - C n-VMR s dC- f Rj . v vjx . n R s dC/2 ( 94 ) 
2 -Xcs J Cj; 


where is the contour of X, the surface of the panel enclosed in the Mach cone, and n is 
the outward normal to the contour lying in the plane of X. 

We now evaluate the integrals in equation (9.4). Since m is constant with respect to the 
variables of integration, the first integral becomes 


- . r e k /* r K(0) 

c f 2 / d e f — ^ 

K S K X Jr, 


0 K _1 *t) \4 2 - r 2 

6 K 


. r 2 = (r?-y ) 2 + (f-z ) 2 


= M S f [x- y x 2 -rK 2 ( 0 )] 

K . 


d 0 


0 K-1 


- 27tCqhx -ju J 
K 


/ 

0 K -1 


0 K 


H K ( 0 )d 0 


(9.5) 
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where Hj^(0) = ^A 2 C q = 1 if the point y, z is inside the panel, and Cq = 0 

for points outside. Similarly, the other volume integral becomes 


/* 0 K r r K (0) 

v.v M J R s d 7 ?df /2 = (v. v M /2) 2 J I /x 2 -r 2 

K /q _ , n 


rdrdfl 


0K-1 0 

0 K 


= (v-'7p/2) 2 [x 3 -H K 3 (0>] 

K , 


d0/3 


0 K-1 


. vp J*R s d7?df/2=jT7. vp 2 ttC 0 x 3 - 2 ^ H K 3 (0)d0 

K „ 


Ok -1 


We have now reduced the potential 0 S to a sum of contour integrals around the area 2 
of the panel contained inside the Mach cone. Since R s (or H) vanishes on the cone, the line 
integrals reduce to integrals over the panel edges contained inside the Mach cone. 

We now evaluate the integration over each edge. For this it is convenient to use 
coordinates oriented normal and tangent to each panel edge as shown in figures 21 and 22. 
It is seen that the angle 6 in terms of these coordinates is given by 


„ i 1 0-y c i -i / 

Vc-u U~ z c , 


from which 


d0=-y c df/r 2 , r 2 = y e 2 + (f-z e ) 2 


(9.6) 
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The derivatives of m in the contour integrals of equation (9.4) become 


VM • n = My 


and 


Rj • v vm • n - (?" z e)^y e z e -ye^y e y e 
Then the last integral of equation (9.4) becomes, for each edge, 

/ R • * • n R s d£ = / [(?-z e ) My e z e - yeMy e y e ] * 

C K 0 


AS 


C K r ~ J!k - 

= My e z e R S 3 / 3 I - yeMy e y e J R 

U •'n 


% 


AS 


since R s 2 - x 2 - (?-z e ) 2 - y e 2 
We also have 


/- 0 K /*% 

J H K ( 0 )dt = -y e y R s d?/r 2 


®K-1 


/ ^K /•*£ /•**■ 

H K 3 (0)d0 = -y e x 2 / R s df/r 2 + y e / R< 

0 K _i ^0 

/. . n\m = Vi yt J' R s d? 






d? 


Utilizing equations (9.5), (9.7), (9.8), and (9.9) in equation (9.4) yields the potential 


0g = Cq [mx + x 3 V. v m/6] + 23 0sk 

k 


(9.7) 


(9.8) 


(9.9) 


( 9 . 10 ) 
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where 



Now since the doublet potential (p^ is given by differentiating 0 S with respect to x 
we obtain 


0d = c O l M + x 2 v. vn/2] + 53 0k 


where 


2 *r0dk = M xy e / df/r 2 R s + v 
*"0 



/* g / 

x3 y e 1 

[ df/r 2 R s - xy e / 

3 * / 0 


d?/R s 


+ 2xy e f R s df/r 2 /6-/u y x /* 

*T> J *T) 

£ /• 2 

d?/R s + XMy e z e R s /2 | + My e y e xy e / df/2R 

0 */ n 


We can express all integrals in terms of two basic functions. Now 
y e f R s df/r 2 = x 2 y e /* df/r 2 R s - y e f 

J 0 J o 0 


d?/R« 


(9.11) 


(9.12) 


(9.13) 


(9.14) 
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We define the basic functions as 


Qi = xy e 



df/r 2 R s ,w n = 



f n df/R s 


The functions Qj and w n are integrated in section A9 of appendix A. 
The doublet potential then can be written 


(9.15) 


0dk = 2^ {mQi - XMy e w 0 + x My e z e R / 2 “ x y e Mz e z e wq/ 2 + x2 Ql(M yy + M zz )/ 2 } (9.16) 


since V-V/i is invariant with rotation of coordinates. 

The source potential is found from equation (9.10) by neglecting the quadratic terms. 
We obtain 


<t>s = CQ a x + £ 0 s k 
k 


where 


<t> sk = oxQi- <7y e [(x 2 ~y e 2 -z e 2 ) 'Vq + 2z e W! - w 2 ] 


The functions W 2 and w j can be expressed in terms of wq by using the recursion 
relations for w n in section A9 of appendix A. This leads to 

0 sk = ctxQj - o ye [ (x 2 -y e 2 )w 0 + (f-z e ) R s ]/2 
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9.2 VELOCITY COMPONENTS FOR THE DOUBLET DISTRIBUTION 
WITHOUT EDGE VORTICES 

When the line integral of the doublet strength is discarded from the integration over the 
panel, the velocity vector is given by the formula 



* 

[n x v 'ju] x y \p drjd? 


where 


(9.17) 


* = l/Rg, R s 2 = (*-x) 2 - (n-y) 2 - (?-z) 2 
v'= a/ag, a/ar?, a+a? 


and 


5= -3/35,3/a 77 , a/a? 


The asterisk denotes finite part of the integral in the sense of Hadamard [11], since 
the integral does not exist in the usual sense. Equation (9.17) can be written in the following 
form after expanding the triple cross product: 




27 rV = / [n • v v'ju - n v'ju . v \ p~\ dr^df 


We assume all of the surface 2 lies in the plane x = 0. Then n = ( 1 , 0, 0), n = jx ( 17 , f ), and 
we have for the x component of velocity 



v f ju • v 1// drjdf 




v'^d^df 


(9.18) 


where now, 


v' = o/ar?, a/a?). 
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For the components lying in the plane of 2 we have, since 9/9£ = - d/9x, 


2ttV= f v-V dr?df 
Z 

Now equation (9.18) may be written in the form 


2n 


/ *§* >• ^ 

v > • v'lpdridt; = - / [v' • ( v'am/ 0 - v ' • v ' ju*/'] dpdf 


Application of the divergence theorem to the first integral yields 


2 tt u = ■ 


f 

Cz 


v/u 


ni/'dfi 


♦ / 


v . v jui//drjd? 


where n now denotes the normal of the contour in the plane of the panel. To express the 
integral in a form which can be integrated without knowing the exact coefficients of ju, 
we expand about the control point projection y, z and obtain 

M 7 ?,?) = M(y,z) + (v - y)M y + (? - z)m z + v yy (i? - y) 2 /2 + M yz (f? - y) (? - z) + mz Z G" - z ) 2 / 2 
For the gradient we have 

v'ju = [My + Myy(n - y) + My Z (? r z)] , [m z + Myz(*J " y) + M ZZ (f “ z)] 
or in vector notation 

v’m =(My +VMy • Rl , Mz + VM Z • Rl) 

where Rl = (V - y, f - z) and v= (9/9y, 9/9z). 

We also have 


(9.19) 


(9.20) 
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yield 


For straight edges, ny and n z are constants and can be moved inside the V operator to 


vju • n = M n +v M n * R 1 


(9.21) 


also 


V . V M = Myy + Mzz = V * V M 


(9.22) 


Equation (9.20) becomes, after substituting equations (9.21) and (9.22) 


27m = - / [/a n + vju n . Ri]i//d£ +v*vn J t//dijd£ 
ci 2 


(9.23) 


Similarly, equation (9.19) yields 


27rv = f +VM y * R G ^ dT?df 


(9.24) 


2 ™ = £- f K + ^-Ri]« 


(9.25) 


where integration is to be performed first. 

Except for the first derivatives with respect to y and z of ju in equations (9.24) and (9.25), 
the integrals are expressed in a form independent of coordinate system. To evaluate the line 
integrals we introduce coordinates related to each edge. First, we perform one integration 
of the area integrals. 



i//d7?df = 


if 

0 K-1 


*0 K 


J rdrdfl/ - r^ 


(9.26) 
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where 6 ^ and are the angles, the kth edge of 2 makes with the Mach cone center y, z. 

Integrating equation (9.26) yields from equation (9.5). 


/* r° K 

J = 27tCqx - 2 J HK(d)d0 

2 K 0K-1 


where Cq = 1 if y, z lies inside 2 and 0 otherwise. 

/ 

Since Rj \jj = -V R s , equation (9.24) becomes 


27TV = 


-bx 


0dr?d? - r Vju y . v R s dr?df 

2 /s 


i 


Since fi does not depend on 77 and f , we have 


27 rv = 

dx 


{"y/ t//dr?df - / v'.(vMyR s )dT?d£'j 


Applying the divergence theorem to the second integral yields 


o 9 

27TV - T— 
dx 


H \pdr}dt; - J* VjUy.'nRgdcj 


dx 


H i/'drjdf K vn-y • nd 2 /R 


From equation (9.27) 


27TV = 


Fy 


27rCo - x 2 / d 0 /Hx( 0 ) 

K 0R-1 

- x I V[iy . nd£/R s 


(9.27) 


(9.28) 
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Similarly, equation (9.23) becomes 


27TU = v • 


VM 


2itCqx - 



H K (0)d0 



+ v/x n • Rj] i//dfi 


The portion of the contour which are arcs of the Mach cone do not contribute to the 
finite part of the integrals since the integrand is infinite on the Mach cone. Restricting the 
integrations over the straight segments allows us to remove the finite part since now the 
integrals exist in the usual mathematical sense. Since R s 2 = x 2 - r 2 , we can write 


y* £ /* £ /• ^ 

J H K (0)d0=-y e J R s d?/r 2 =-y e x 2 J ^+y e J ^ 


^K-l 


0 ‘ " s 

= - xQi + y e wo 


( 9 . 29 ) 


Since Mn — My e 


/ C 

[- yeMy e y e + (f-z e ) My e z e ]dr/v/x 2 -y e 2 -(r-z) 2 




--yeMy e y e wo-R s My e z e 


Finally, the velocity components can be written 


u = xv*vpCo + 2 u k 

k 

v = M y Co + S v k 

k 

w = M z Co + 2 w k 
k 


( 9 . 30 ) 
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where 


u k = [- My e wo - yeMz e z e wo + My e z e R s + xQjv. vm]/2tt 
v k = [^yQl - XMyy e W 0 ] /2 tt 
W k = [m z Qi ~ x Mzy e w 0 ]/2tt 


9.3 VELOCITY COMPONENTS FOR SOURCE DISTRIBUTION 

The velocity components from the source distribution on a superinclined panel are 
obtained by integrating the gradient of equation (9.1) or 


2tt V = 



o(-o£)\p dfjdf 


where v = 9/9x , 9/9y , 9/dz. _ _ 

It is convenient to separate the velocity vector into components u and V pj where V p is 
the velocity vector in the plane of 2. We easily see that 


ai//dr?df = </> d 


where d> d is the potential for the doublet after dropping the quadratic terms. Now Vp may 
be written 



2ttV, 


v <^di7df 


= - J v' (o\p) drjdf + J v' at//d??df 


where v’ =(9/9 tj , 9/9 f). 

Applying the divergence theorem to the first term yields 


' v p- / 


2tt V n = - J na (rj , f) * dC + va / ^df 

Cs 2 


"I 


(9.31) 


(9.32) 


(9.33) 


(9.34) 


(9.35) 
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Here we have used the property that o is linear. 


Now the curved parts of the contour make no contribution to the finite part. Restricting 
the contour integration to the straight segments allows us to eliminate the finite part since 
the remaining integrals exist. In terms of the coordinates y e , z e associated with the edge , 
the expansion of the source distribution o takes the form 


a(V , ?) = ff(y e , z e ) + Oy e (0-y e ) + o Zq (f-z e ) (9.36) 


Then the first integral in equation (9.35) becomes, for each straight segment 

r g A f £ (f-zp’) dC 

- ( ff - ff y e ye) "k / Vx2 - y e 2 - (J-z e ) 2 - o Zq n k J V X 2 - y 0 2 - (f-z e ) 2 

*'o 0 

= [- ( a - a y e ye) wo + a Ze R s ] n k (9.37) 

where functions wq and R s are evaluated at the endpoints of the side, or its intersection with 
the Mach cone. 

The second integral in equation (9.35) has already been evaluated in equation (9.27) 
and (9.29). The u component is found from equations (9.12) and (9.16) after dropping the 
quadratic terms. Hence, the velocity components take the form 


u = Cqcj + Z u^ 

V p = C 0 V(7 + S V pk 

where 

u k = [ CT Ql “ xa y e wo]/ 27 r 

Vp k = n k [(y e Oy e - o) w 0 + a Ze R s ] /2 tt + v a (xQi - y e w 0 ) /2 tt 
and n k is the outward normal for the kth edge in the plane of the panel. 


(9.38) 


(9.39) 
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10.0 TRANSFORMATION OF VELOCITY COMPONENTS FROM 
LOCAL PANEL TO THE GLOBAL COORDINATE SYSTEM 


The influence coefficients for a given panel have been obtained in terms of local scaled 
panel coordinates. To find the contribution to the induced velocity on a given panel from 
all panels of the curved surface to be analyzed, we use a transformation from each panel 
coordinate system to the global system. For this purpose, we consider three basic coordinate 
systems: 

1. Reference coordinates (x r , y r , z r )chosen for convenience in defining the body surface 
grid on which the panels are constructed. 

2. Compressible coordinates (x 0 y c , z c ) in which the velocity of the free stream is aligned 
with the positive x c axis and, of course, 

3. The local scaled panel coordinate system ^Xp, yp, Zp) 

10.1 TRANSFORMATION FROM THE REFERENCE COORDINATE SYSTEM 
TO THE COMPRESSIBLE SYSTEM 

It is frequently convenient on a configuration to use coordinates in which the x axis 
is not aligned in the free stream direction. We then must establish a transformation from 
this coordinate system to the system in which the compressibility scaling is applied. From 
the reference coordinate system x r , y r , z r , we apply a rotation of the velocity vector through 
an angle of yaw j3 c followed by a rotation through an angle of attack a c (see figure 23). 

Here a c and j3 c are not necessarily small. The first transformation in figure 23 is a rotation about 
the y r axis, Thus 


x 0 = x r c os/3 c + y r sin/3 c 
Y 0 = -x r sinj3 c + y r cos i/3 c 


z 0 z r 


( 10 . 1 ) 


We then apply a rotation a c about the yg axis. This leads to the following transforma- 
tion to compressible coordinates x c , y c , z c : 


x c = xq cosa c - zq sina c 

y c = yo 

z c = XQ sina c + zq cos a c (10.2) 


70 



Combining the results of the transformations yields 


x c = x r coso: c cos/3 c + y r cosor c sinj3 c + z r sina c 

y c = - x r sin/3 c + y r cos/3 c 

z c = " x r sina c cos/3 c - y r sina c sinj3 c + z r coso: c 


On defining the transformation matrix | A r J equation (10.3) takes the form 


( x c> y c > z c) {^r} ( Yr 


where 


Kl 


cosa c cosj3 c cosa c sin/3 c sina c 


sinj3 c 


cosfL 


0 


, - sino: c cosj3 c - sina c sin/3 c coso: c , 


When local panel coordinates Xp, yp, Zp are related to compressible coordinates by 

/ X C N 

( X P> y p . Zp) = {a c [[ y c 


we have 


( x P’ y P 5 Zp) - |A C | |A r | j y r 


(10.3) 


(10.4) 


(10.5) 


( 10 . 6 ) 
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where { A c } is the matrix of transformation from the compressible coordinates to panel 
coordinates. Solving the set of equations (10.3) for x r? y r , z r in terms of x c , y c , z c yields 

x r = x c cosa c cos Pc “ Yc ” z c si na c cos Pc 
y r = x c coso: c sinj3 c + y c cos/3 c - z c sina c sinj3 c 

z r = x c si ncv c + z c cosa c (10.7) 

and we see that 


/cosa c cosp c - sin] 3 C - sina c cosj3 c ^ 

■{A r 'l} = [ cosa c sinj3 c cos]3 c - sina c sin/3 c 


smo: r 


and 


cosj3 c 


( x p y r’ z r) {^r y c ] H 



( 10 . 8 ) 


(10.9) 


The normal in the compressibility coordinate system is defined in terms of the reference 
normal by 


{ A r [ fi=(\0).fi,A r <2).a, A r (3) . n) 


( 10 . 10 ) 


where A r (i) is the ith row of the matrix | A r | . The conormal in the compressibility coordinate 
system is 


(±B2 Ar (l).a, A r (2).ft, A r (3).ft) 

where the plus sign holds for subsonic flow and minus sign for supersonic. 


( 10 . 11 ) 
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If we define the matrix ) B ) by 



( 10 . 12 ) 


then the conormal in the compressible system may be expressed as a matrix product 

)B( )A r )fi 

In the reference coordinate system the conormal vector n c is given by 

n c T = {A r -1} {B} {A r f ft 


(10.13) 


(10.14) 


where n c ^ is a row vector. 

10.2 TRANSFORMATION FROM COMPRESSIBLE COORDINATES TO PANEL 
COORDINATE SYSTEM WITHOUT THE COMPRESSIBILITY FACTOR B 

IN THE PANEL SYSTEM 

For incompressible flow, the differential equation and distance are invariant under 
rotations about any axis, while for supersonic linearized theory, the differential equation and 
the hyperbolic distance are invariant under rotations only about the free stream direction 
(x c axis) and under any oblique or Lorentz transformation in other planes. It is convenient 
to introduce at this point the transformation from scaled coordinates to compressible 
coordinates; namely 


(x s , y s> z s ) = (x c , By c , Bz c ) 


(10.15) 


We now consider a rotation through an angle 6 j about the streamwise scaled com- 
pressible coordinate axis. This becomes 

xi = x s 

y l =y s cos0j -z s sin0 j 

zi = y s sin0 j + z s c°s0 1 (10.16) 
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Introducing compressible coordinates into equation (10.16) leads to 

xi =x c 

y i = By c cosfl j - Bz c sind j 

z = By c sind j + Bz c cos(9 j (10.17) 


We chose the angle d j so that the normal lies in the y j = 0 plane (see figure 24). Then 

tan 0 1 = y c /z c = ny c /n Z( , (10.18) 

from which sin d j = n-yj y/n^- + n Zc 2 and cosd i = n Z J >/ny c 2 + n Zc 2 . 

The normal will be calculated generally in the reference coordinate system. Hence, we have, 
for the normal in the compressible system, from equation (10.10) 

n x c = A r O ) . n , n Yc = A r (2) . n, n z ^ = A r ^) • n (1 0. 1 9) 

where A f W denotes the ith row of the matrix j ^ [ 

For subsonic flow, we rotate the xj, yj, z\ coordinates system about the y\ axis and 
the transformation becomes 

x 2 = x i cos^2 " Z 1 sin^2 

y 2 = yi (10.20) 

z 2 = X 1 sin^ 2 + Z 1 c °s^2 

We chose 6 2 to make the panel lie in the Z2 = 0 plane; then from figure 25 


tan #2 - - z 1 /x 1 = n x 




( 10 . 21 ) 


from which 


sin02 = n x an< ^ cos ^2 = V^ 


y s 2 + nz s 


( 10 . 22 ) 


Since the normal in the scaled coordinate system is related to the compressible system by 

= B n x / N /B2n x 2 + n y 2 + n z 2 
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n y s n yJ + n y c 2 + nz c 2 


n z s = n z c /\/ fi2n x c 2 + n y c 2 + n z c 2 

(10.23) 

then 

sin0 2 = Bn X( ,/ v/B 2 n Xc 2 + n y 2 + n z 2 


cosd 2 = \/ny c 2 + n Zc 2 / v/B 2 n X(; 2 + hy c 2 + n Z( , 2 


Since the vector B^n v m, , n 7 is the conormal, we have 
A c’ ^c ^c 


sin# 2 = Bn Xc> /n • n c 

(10.24) 

cos0 2 = \/ n y c 2 + n y c 2 /\Z" . n c 

(10.25) 

For supersonic flow we apply an oblique transformation in the x\ , 
general form 

z j plane. This has the 

x 2 = X 1 /v^l - nii 2 + mjzj/v/l - mj 2 


y2 = vi 


Z2 = mixi/ v/l - mi^ + zi/^/l - 

(10.26) 

For subinclined panels, we chose the coordinate Z 2 = 0 as the panel, then we obtain from 
figure 25 

mi = - zi/xi = tan# 2 

(10.27) 
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Now 


tan 02 = 



= m l 


Since 


n Y ~ Bn Y 
A s A c 

n y s ~ n y c 

n z s ~ n y c 


we have 




Bn x c 

V 7 ' - n x c 2 


= m l 


from which 




/ y/l ni]2 Bn x< J v^y c ^ ~^ n z r ^ ~ ®^ n Xr ^ 


= Bn x / v/ft- n c 


and 


l/v/l — m j 2 = v / ny c 2 + n Zr 2/ v /^T^ 


where n c is the conormal given by 


n c =(-Bn Xc 2, n^, n Zc ) 

and 

" = ( n x C ’ n y C ’ n ^c) 


( 10 . 28 ) 


( 10 . 29 ) 


( 10 . 30 ) 
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Combining the two transformations yields the following tranformation from compressible 
coordinates, x c , y c , z c to x p , y p , z p = X 2 , Z 2 J coordinates for supersonic flow: 

x p = x c/ \A + m l B (y c sin0 j + z c cosfl j) / >/l - mj 2 

y p = B (y c cosfl \ - z c sinfl j ) 

z p = m i x c / - mj 2 + B (y c sin # \ + z c cos# / ^/l - mj~2 (10.31) 

The determinant of the matrix of equation (10.31) is equal to B^. Solving for x c , y c , z Q 
in terms of x p , y p , z p yields 

x c = x p/v 1 -mi 2 -m 1 Zp/ v 'l - rri]2 

y c = - Xpm] sin0 i/B-y/l - mj 2 + y p cos0 j/B + z p sii^i/B-y/l - mj 2 

z c = -x p mi cos0 j/B^/l -mj 2 -y p sin0j/B + z p cos0 i/B^/l^mT 2 (10.32) 

For superinclined panels we choose mj so that X 2 = 0 lies on the panel, then from 
equation (10.26) we get 

m l = - x \/ z \ ~ cos#2 

m l = \A ~ n x c ^ /Bn Xc = y/\ -^Aj-0) • n)2 ( BA r O) • n c (10.33) 

from which 

mi/v/l -m 1 2= \/ny c 2 + n Zc 2 / N /-n • n c 

and 

l/v/l -m 1 2 = Bn Xc / v /-n. rT c 

Since I m j I < 1 we see that 

!/ B < n Xc /v^ = tan 0 2c 
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and the angle 02c between the panel and free stream must be greater than the Mach angle. 

The matrix and its inverse have the same form for both superinclined and subinclined panels, 
but differ in the calculation of the parameter m j. 

For the sake of convenience in the combined program for both subsonic and supersonic 
flow, we make the transformation for subsonic flow and for the subinclined and superinclined 
panels take the same form. Hence, we must rotate the coordinates of the superinclined panels 
so that Xp, y p are coordinates in the panel. Then the new panel coordinates are 


x p = z 2 

y p = -V2 

z p = x 2 (10.34) 


The variables X 2 , Y2> z 2 correspond to x, y, z in the derivation of the aerodynamic influence 
coefficients. Note that for superinclined panels the value of m j is the reciprocal of the value 
for subinclined panels. Let m 2 be the value for superinclined panels; then substituting equa- 
tion (10.34) and replacing mj by m 2 = 1/mj into equation (10.29) yields 


X P = X 1 / v/in 1 2 _ 1 +mizi/ v /m 1 2 - 1 
y p = -Yl =yi Sign (l -mj2) 

z p = mixi/ v /m 1 2- 1 +z 1 / v / m 1 2- 1 ( 10 . 35 ) 


with m j defined by equation ( 1 0.30). When m j ^ - 1 is replaced by | 1 - m j ^ | then equation 
(10.35) holds for both subinclined and superinclined panels. Since for subsonic flow 


B 2 n Xc 2 + n Yc 2 + n Zc 2 = n . n c 


(10.36) 


we find that the transformation from the reference coordinate system to the panel coor- 
dinate system in both subsonic and supersonic flows can be written in the form 


fx r — 


x 0> 


( x p> y p > z p) = {a 2 } {A!} |A r } [ y r -yo 

k z r - z 0 4 


(10.37) 
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where { A r | is given in equation (10.5), xq, yq, zq is the center of the panel and the origin 
of the panel coordinate system. 


and 



0 

Bn zJ7 


En yjy Bn z c h 



(10.38) 



s Bn Xc /V|n • n c | 


sin (n • n^) 


Bn x r /v / |n • n c | 0 



(10.39) 


where s = -1 if the flow is supersonic and s = 1 if the flow is subsonic; 

and B-yH ~M 2 I. 




7=V n “y c + n " 


(10.40) 


The points defining the panel are required in the average plane panel coordinates. 
Similarly the corner points defining the triangular subpanels are also required in the local 
subpanel coordinate system. All these coordinates are computed using equation (10.37) as 
defined for the particular panel or subpanel. The projection of the curved panel on the aver- 
age plane is found by setting Zp = 0. 

For supersonic flow, the quantity n . n c can be seen to be a test for sub- or superinclined 
panels. If n * A n c J> 0, the panel is subinclined, if n • n Q = 0, the panel is inclined at the Mach 
angle; and if n • n c < 0, the panel is superinclined. Panels inclined at the Mach angle are not 
admissible since the transformation to compressible coordinates breaks down. This can be 
seen by setting ni| = 1 in equation (10.26). 

10.3 CORRECTION TO PANEL ELEMENT OF AREA IN SUPERSONIC FLOW 
DUE TO THE TRANSFORMATION FROM GLOBAL TO 
PANEL COORDINATES 


The influence coefficients are computed in the coordinate system of the panel inducing 
the flow. Hence the conormal derivative in the compressible direction and the element of 
area ds must be transformed. The factor from page 383 of Courant [16] , Vol. 2, relating the 
surface panel element of area dXpdyp to ds is 


ds 


= n > 


d(y C * z c) 9 ( Z C. x c) d( x c> y C ) 

d ( x p’ y p) nyc a ( x p> y p) nzc 0 ( x p’ y p)_ 


dxpdy p 


(10.41) 
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Evaluating the Jacobians from equation (10.32) yields 


ds = [n Xc m 1 /B+ (n^sin^ + n Zc cos0 } ) ] dx p dy p /B yf\ -mj 2 
Substituting for m j from equation (10.30) and 6 j from equation (10.18) finally gives us 


ds = dXpdyp/B^/l - M 2 n Xc 2 


(10.42) 


for subinclined panels. For superinclined panels the area element on the panel is in the y p , 


Zp coordinates and surface elements are related by 


ds = 


3(y c > z c) . 9 ( z c> x c) , 9 ( x c> y c ) 

" X C 9 (y p , z p ) n yc 9(y p , z p ) nz c 9(y p , z p ) 


dy D dz 


p^p 


(10.43) 


Here we used the transformation in equations (10.31) and (10.32) instead of equation (10.34) 
through (10.39). Evaluating the Jacobians using equation (10.32) yields 

ds = [n Xc /B + m l ( n y c sinfl \ + n z c cosfl j) j dypdZp/B^/l - mj 2 
Finally, substituting equation (10.33) for m j and equation (10.18) for gives us 

ds = dy p dZp/B - 1 ( 1 0.44) 


for superinclined panels. 

The integrals involving the doublet do not need the area factor if the conormal deri- 
vative is replaced by the derivative in the panel coordinate z p for subinclined panels and x p 
for superinclined panels. For subinclined panels, the conormal derivative takes the form 


0$/3 n c = - B^n Xc 00/0x c + ny c 00/0 y c + n Zc 00/0 z c 

9<i/9n c = B2n Xc 9x p /9x c + n yc 9x p /9y c + n Zc 9x p /9z c ) 90/9x p 

+ (- B 2 n Xc 9y p /9x c + n yc 9y p /9y c + n Zc 9y p /9z c ) 9<£/9y p 
+ (- B 2 n X( , 9z p /9x c + n yc 9z p /9y c + n Zc 9z p /9z c ^90/9z p 
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The coefficients of 


30 

3x r 



are found to equal zero when the derivatives inside the 


parentheses are evaluated from equation (10.31) with 8 j and mj determined by equations 
(10.18) and (10.30), respectively. Evaluating the coefficient of 30/3z p yields 


80 

9n c 


= -B> 2 „x c 2 -‘f^ 


(10.46) 


and we see that 


80 , _ 80 , , 
9n c ds “ 9zp dx P dy P 


(10.47) 


When the same procedure is applied to the superinclined panels we obtain 

!*- = B^l -M2 „ x 2 |4. (10.48) 

8n c v x c 3Zp 

and 

^0 . 30 „ J 

9^ ds = “9^ dy P dz P ( 10 - 49 > 

Hence the area correction term need not be applied to the doublet integrals. 

10.4 TRANSFORMATION OF VELOCITY COMPONENTS 
TO COMPONENTS IN THE REFERENCE SYSTEM 

The velocity components are computed in the coordinate system of the panel inducing 
the flow. In terms of the compressibility direction, the velocity components are found from 
the chain rule; i.e., 


80 _ 80 d x p 00 3 y p 00 8zp 

3x c 0Xp 0x c 8y p 3x c 0Zp 8x c 


(10.50) 


00 00 

with similar relations for — — and -r— . 

9y r 9z r 


90 90 

^ ^ > and 

9x p 9y p 


30 

3z„ 


yields 


Using equation (10.31) to evaluate the coefficients of 
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30 

_ 00 

... 1 , 

00 

mi 



dx c 

3x r 

\/l - mj 2 

3 z p 

0 -mj 

7 


3 0 

00 

m B si nO 

00 


00 

B sin0 

3y c " 

0X p 

x/l - mj^ 

9y p 

Lj cost / 1 

0Zp 

1 y/\ - m^ 

30 

_ 00 

m i B cos 6 ] 

30 


00 

B cos0 \ 

dz c 

“0X p 

y/l-nfl 

0y c 

d sine / 1 + 

3 z p 

yj\ - mj 2 


We see that this may be expressed in the form 

v c = |aT}v p ( 10 . 52 ) 

where {aT} is the transpose of {A} , the matrix of equation (10.31). 

We now resolve the vector V c into components along the coordinate axes of the reference 
coordinate system. Thus 


v- )A r ->|V c HA r T|V c = )A r T( | A T(v p 

or 

V = (l A i iM) 1 V p (10.53) 

This formula is valid for both subsonic and supersonic flow. 
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1 1 .0 FAR FIELD EXPANSIONS FOR THE 
AERODYNAMIC INFLUENCE COEFFICIENTS 


11.1 SOURCE AERODYNAMIC INFLUENCE COEFFICIENTS 

For panels lying wholly within and sufficiently far from the Mach cone emanating up- 
stream of the control point, far field expansions of the aerodynamic influence coefficients 
are feasible. These expansions are algebraic and are much more economical in computing 
time. Except for a few signs, the subsonic and supersonic far field expansions are identical 
in form and a single derivation of the far field aerodynamic influence coefficients is possible . 

It is convenient to represent the far field integrals in the compressible coordinate system 
rather than in the panel system. We shall translate the origin of coordinates to the center of the 
panel inducing the flow to simplify the analysis. Let P denote the control point and Q be 
the location of the singularity point on the panel and the coordinates of integration. The dis- 
tance, or hyperbolic distance, is given by the scalar product 

Rb = [p-Q.p-Q ] 1/2 (ii-D 

where the scalar product [A, B] is defined by 

[A, B] = ajCjjbj 

and Cjj is the matrix 

/. 0 

{c} =Cjj= I 0 sB2 

\0 0 


and B =>/jl - M^| , s = 1 for subsonic flow and -1 for supersonic. 
* 

To utilize subscript notation effectively, we let 


( 11 . 2 ) 



P= (xj, X 2 , x 3 ) 

and 

Q = (y1.y2.y3) 
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The expansion of 1/Rg then becomes 


_U_L_ y ._L 
R b R c Yl 


/I VI 9 2 
(rc/ 2 y » y j axj axj 



1 

Rb 


1 . yj c i i x i 

R c R c 3 



c ij , 3xjXj\ 

Rc 3 Rc 3 / 


-•^yiyjyk ( 3 /R c 5 ) ( c ik x j +Cjk~i+CyX k )-15x i xjX k /R c 7 


1/2 

where R c = [p, p] is the distance for subsonic flow, or hyperbolic distance for super- 
sonic flow, from the panel center to the control point and Xj = XjCy. Using the fact 
that yjX^ = y^X|, we obtain 


^ = GW) + y iGi U) + G U <2) +-Lyfty k G 0k W 

where 

G(0) = 1/R C 
G/D = Xj/R c 3 

G ij (2) =3x i xj/R c 5 -c i j- 1 /R c 3 

G ijk (4) = 15x i x j x k/ R c 7 _ ( 3 / R c 5 ) ( c ik" lx j + c jk" lx i + c ij _lx k) 


To compute the velocity vector we obtain 


| v p (1/ R b)| k a x k ( R c) y ‘ axj dxR (r c ) + 2 Yiy i axj 0Xj 0XR (r c ) 


{- c _1 vp (1 /Rb) } k - G k ( 1 } + yi°ik (2) + G ijk (4) + • • • 


(11.4) 


(11.5) 


( 11 . 6 ) 


(11.7) 



The source potential is given by substituting equation (1 1.5) into 


0 = -( 1/co) 


/ 


ads/Rg 


where co = 2ir for supersonic flow and 4ir for subsonic. We obtain 


0 S = - G(O)Ed) - GjCDe^) - G ij (2)E ij (4) 

where 

Ed) = (1/co) / ads 
Ej(2) =(l/a>) J' ay ids 

Ejj( 4 ) = (l/a>) / ayjyjds 


Similarly, the source velocity components are given by 

V.-i/ v p (l/RB)ds 

ic-n v s =i / -[-{c-n v P (i/R B )]ds 

{C- 1 } V s = G k (l)E(l) + G ik (2)Ei(2) + Gy k (4)Ey(4) 


The integration of the E integrals are performed in panel coordinates, 
of the y variables to panel coordinates is affected by 

Y= {C} Y= {C} {a -1 } (£, tj» 0) 


and the element of area is 


ds = J d£ drj 


( 11 . 8 ) 


(11.9) 


( 11 . 10 ) 

The transformation 
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where J is the Jacobian of the transformation from compressible to panel coordinates whose 
matrix is A. Since o = o Q + cr x £ + o^r], the Matrices E are products of the matrices relating 

°o> a x s a v *° va l ues °f ^e source at panel centers of neighboring panels, of the matrices 

{c} , { A“^} , and the moments 


Cmn’ff ( m+] 


r? n+1 d£dr? 


( 11 . 11 ) 


where integration is over the panel. The development of the formulae is fairly straight for- 
ward but tedious and is omitted. 

11.2 DOUBLET AERODYNAMIC INFLUENCE COEFFICIENTS 

The doublet potential is given by 


<t>=lo U * V P 0/ R B) ds 

~ ‘oJ J" V n c {g} v p (l/ R B)J ds 

Since n c = n{B} , then n c |c| = jiT{b} jcj- = sB^n^ . Here 





and n^ is a column vector, 
then 


where 


0 = " G k ( 1 >W k ( 1 ) - G ki (2)w ki (2) - G kij (4)w kij (4) 

W k (D =(sB2/cj) i J /xn k ds 
W k j(2) =(sB2/co)^* pn k yjds 
w kij (4) = (sB 2 /co) J pn k yjyj ds 


( 11 . 12 ) 


(11.13) 
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The evaluation of the W matrices is carried out in a manner similar to the E matrices. 


The velocity components for the far field are evaluated with the line integrals of the 
doublet strength neglected. One can easily verify that in the second term of equation (3.51) 

B2(R 0 -R)/R B 3 = _ {b} v p (1/R B ) 

= - {bHc} [- lc-q v p (i/r b )] 

= sB 2[-{ C -1} v p (l/R B )] (11.14) 


Substituting the equation (11.14) into the second term of equation (3.5 1) yields, with the aid 
of equation (1 1.7), 

V d = (nXvp) X [- {c -1 } Vp(l/R B )]ds 

= (n Xv M ) X [g( 1 ) + y iGi (2) + \ y^Gy^)] ds 

V d = Gd) xF(l) + Gj(2) XFj(2) + Gjj (4) XFy(4) 


(11.15) 

(11.16) 


The sign of equation (3.51) has been changed since we are using a normal n into the fluid. 
Here 


F(l) = (sB 2 /co) ^(nXvju)ds 

^i(2) = (sB^/cj) J* (nXVM)yids 

Fy(4) = (sB2/co) J (n XVM)~iyjds (11.17) 


Again, the integration is over the panel. The F matrices are again products of several matrices, 
relating transformation to panel coordinates, relating the coefficients of p to surrounding cen- 
ter panel values of doublet strength, and of the moment integrals C mn in Equation (11.11). 
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12.0 DISCUSSION OF RESULTS 


12.1 SOURCE PANEL METHOD APPLIED TO CONES AND 
BODIES OF REVOLUTION 

The source panel method was first applied to cones and slender bodies of revolution. 

For a free stream Mach number ofy/2, figure 26 shows the pressure coefficient versus cone 
half-angle for zero angle of attack. The pressure coefficient agrees well with the exact solu- 
tion and begins to depart from the theory only for cone half-angles approaching the shock 
detachment angle. Figure 27 gives a comparison of the exact solution from references [17] 
through [20] for the flow over 10° and 1 5° half-angle cones over a range of Mach numbers 
with the present method, and with the method modified to apply tangential velocity boundary 
conditions instead of tangential mass flux boundary conditions. Solutions with tangential 
velocity boundary conditions generally underpredict the pressure magnitude. The agreement 
with exact theory is very good for the 10° cone up to a free stream Mach number of 3. For 
the 15° cone, agreement with exact theory is excellent up to a free stream Mach number of 
1.6, but drops rapidly for higher Mach numbers. The rapid drop in pressure occurs for near 
stagnation regions of the flow when mass flux boundary conditions are used. A modification 
of the present calculation to give better results in stagnation flow is presented in the section 

12.2 on elliptic cones. 

For a 10° half-angle cone at 5° angle of attack the source panel method with tangential 
mass flux boundary conditions using 36 panels around 180° of the cone axis is compared with 
the solution from tangential velocity boundary conditions and with exact theory of references 
[17] to [20] in Figure 28. The velocity boundary condition underpredicts the pressure espe- 
cially on the high pressure side of the cone. 

Comparison of the pressure coefficient on an axially symmetric spindle computed by the 
linear source panel method with results from the exact method of characteristics and with a 
line source distribution solution is presented in figure 29 and the three methods are in good 
agreement. Here mass flux and tangential velocity boundary conditions would be expected 
to be in closer agreement, since the body slopes are small. Also shown on the graph are results 
from a constant source panel method using the same panel distribution. The present method 
with linearly varying source distributions is seen to yield better results. 

To test the stability of the source numerical method a spindle was paneled in a random 
manner as shown in figure 30. The resulting values of the pressure distribution at the control 
points of the panel is indicated by the dots in figure 31. Considering the extreme variation in 
panel size, we see that the agreement with the exact calculations is very good and that the 
numerical method is remarkably insensitive to panel shape and size. 

12.2 SOURCE PANEL METHOD APPLIED TO ELLIPTIC CONES 

Flows around non-axially-symmetric three-dimensional bodies in the form of elliptic 
cones were calculated using the source panel method. The second order theory of Van Dyke 
[21 ] was programmed and the pressure coefficient computed at the control points to compare 
with the source panel calculated values. Figure 32 shows the two pressure distributions on an 
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elliptic cone with a maximum cone half-angle of 14° and a ratio of minor to major axis of 
0.532 at a Mach number of y/2. Agreement is seen to be very good with the Van Dyke method 
yielding slightly higher values. Similar calculations are shown in Figure 33 for a maximum 
cone half-angle of 15°, a ratio of minor to major axis of 0.3 at a Mach number ofy/T. The 
agreement is excellent. For this case the source panel method gives slightly higher values on 
the flat side of the cone than the Van Dyke second order theory. Similar results are shown 
in figure 34 for a cone with a maximum half-angle of 18.4° and a ratio of minor to major axis 
of 1 /3 at a Mach number of 1.81. Figure 35 shows a graph of the pressure coefficient from a 
cone with a maximum half-angle of 30° and a ratio of minor to major axis of 0.2. Here, the 
panel method yields considerably lower pressure than the second order theory. 

For the cone of 30° maximum half-angle in figure 35 elliptic coordinates were used to 
obtain a finer paneling in the higher pressure region. For this paneling, a sharp dip in the pres- 
sure occurred in the region of maximum cone half-angle_(see figure 36). This dip can be 
remedied by choosing the perturbation velocity vector V to be proportional to perturbation 
mass flux vector W by the relation 


V = W/p, with p = 1 - M^ (j) x + . . . for 0 X < 0. 


This equals the correct velocity to the first order in the perturbation potential 0 and when 
-0 X gets large near stagnation points it gives better pressure results. The pressure distribution 
corrected by this formula using the isentropic pressure relation is shown by the triangles in 
figure 36. The pressure distribution is considerably improved, although the comparison with 
second order theory is not as good as for the more slender cones. 

Figure 37 compares the calculated pressures obtained by using mass flux boundary condi- 
tions with those obtained by prescribing zero normal velocity boundary conditions and with 
experiments. For this example, velocity boundary conditions yield better agreement with 
experiments (ref [22] ). The linear solution with mass flux boundary conditions overpredict 
the pressure particularly on the flat side of the cone. Furthermore, the Van Dyke [21 ] second 
order theory tends to yield higher pressure as in the previous cases, than linearized theory 
with mass flux boundary conditions and, hence, also over preducts the actual pressures on 
elliptic cones. 

12.3 APPLICATION OF SOURCE AND DOUBLET PANEL METHODS 
TO PLANAR WINGS WITH LINEARIZED BOUNDARY CONDITIONS 

Source and doublet panels distributed on a wing planform were used to represent thick- 
ness and camber with linearized boundary conditions. Figure 38 shows the pressure coefficients 
for a yawed flat plate delta wing at angle of attack computed by the present method compared 
with the exact linearized theory solution in Jones and Cohen [23] . Agreement is seen to be 
good. However, the discontinuity in velocity gradient on the Mach line through the apex can 
be represented even better by a paneling which uses the Mach line emanating from the vortex 
as a network boundary. For this case (not shown in the figure) the pressure was exactly con- 
stant in the region between the supersonic leading edge and the Mach line through the comer 
as the theory predicts. 
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The planar doublet panel method was tried on a wing of zero thickness with a parabolic 
arc camber profile described in figure 39. The pressure is in excellent agreement with exact 
linarized theory solution. Breaking the wing into networks bounded by Mach lines emanating 
from discontinuities in the wing leading edge allows for better representation of the infinte 
velocity gradient occurring at such Mach lines. 

Using the same parabolic arc profile wing, the targential free stream perturbation velocity 
calculated from the theoretical solution was applied as design boundary conditions to a section 
of the wing between the supersonic leading edge and a radial line passing through the comer 
with the subsonic leading edge as seen in figure 40. Analysis boundary conditions were applied 
to the remaining portion of the adjacent to the subsonic leading edge. The computed down- 
wash from the panel method is compared with the actual wing slopes in figures 40 and 41 along 
the two lines numbered N = 1 and 2 in the sketch. Figure 40 shows the greatest departure of 
the downwash from the actual wing slopes, and the match with actual slopes improves as one 
moves closer to the supersonic leading edge. It was found that the design method was unstable 
when the subsonic leading edge was included in the design network. However, these results 
represent only an initial attempt, and it is anticipated that these problems can be overcome. 

With the same planform as in figure 39 and a symmetrical profile described by the same 
parabolic arc, the pressure distribution was obtained by the planar source panel method. Com- 
parison with the exact solution is presented in figure 42. The calculated results are in excellent 
agreement with the theory. Using the theoretical pressure distribution in the form of the 
tangential free stream perturbation velocity component as design-type boundary conditions 
in this case, yielded a downwash distribution which was undestinguishable from the actual 
wing slopes. 

While using planar panels with linearized boundary conditions yields good comparison 
with exact linear theory solutions, a better comparison with experiment and with more exact 
supersonic solutions is obtained by paneling the wing surface and applying the exact boundary 
conditions of vanishing normal mass flux. Figure 43 compares the pressure distribution on a 
symmetric parabolic arc two dimensional airfoil by the present source panel method adapted 
for two dimensional flow with exact linear theory, the second order expansion (ref. [24] ), 
and also with experiment. The present method with the pressure computed by the isentropic 
relation is in close agreement with second order theory on the upper surface and with experi- 
ment. The agreement on the lower surface is not quite so good but is still considerably better 
than linear theory. 

12.4 APPLICATION OF COMBINED SOURCE AND DOUBLET PANELS 
WITH POTENTIAL BOUNDARY CONDITIONS TO 
THREE DIMENSIONAL CONFIGURATIONS 

For configurations of irregular shape such as inlets, nacelles, and bodies, discontinuities 
in slope produce disturbances in the interior which reflect from the singularity sheet represent- 
ing the configuration surface. These waves on repeated reflection may build up in strength 
and affect the surface singularity distribution which will in turn, introduce perturbations in 
the exterior pressure. These perturbations are spurious in that they are not properly .related 
to physical conditions on the exterior surface where they occur. By using combined source 
and doublet paneling with both interior and exterior boundary conditions, the interior flow 
perturbations may be eliminated. Panel spacing requirements then depend only on the exterior 
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geometry as does the true physical flow, without further spacing requirement arising from 
interior reflected waves. 

For a closed body, this may be accomplished by setting the source distribution a equal to 

A 

o - - U • n 


where U is the free stream velocity vector and n is the outward normal to the panel. The 
doublet distribution is then found by setting the perturbation velocity potential 0 on the 
interior surface equal to zero at the panel control points. 

A simple axially symmetric inlet was first analyzed by the source panel method. The 
pressure distribution along with the drawing of the inlet with its panel distribution is shown 
in figure 44. The result indicates that internal waves from the lip are influencing the pressure 
distribution on the outside of the nacelle. To eliminate the internal waves, a combined source 
and doublet distribution was applied to each panel. For this particular calculation, the 0 = 0 
interior boundary condition was replaced by the equivalent boundary condition of the 
vanishing of the normal component of the mass flux to the exterior surface. (Although the 
inlet is not closed, the 0 = 0 interior boundary condition would work successfully for the iso- 
lated nacelle since 0 is identically zero on the inlet opening in the absence of an incident 
perturbation flow.) The resulting pressure distribution, shown in figure 45, is much smoother 
and in better agreement with the Lighthill [25] solution. 

The combined source-doublet paneling method with exact boundary conditions was 
also applied to the Carlson [26] , [27] wing 2T. The paneling is shown in figure 46. The 
pressure distribution on both upper and lower surfaces by the exact isentropic formula is 
plotted in figures 47 and 48 for four of the spanwise locations in figure 47. Except in the tip 
trailing edge region the pressure distribution is smooth and compares favorably with experiment. 
The smoothness of the pressure coefficient indicates that with the combined source and doub- 
let paneling there is little buildup of internal waves. (As we shall see later, the tip region pres- 
sures are improved with the nine parameter spline.) 

12.5 APPLICATION OF THE CONTINUOUS DOUBLET SPLINE 
TO THREE DIMENSIONAL CONFIGURATIONS 

The doublet spline used on each of the preceding examples was developed earlier for 
subsonic flow and is described in reference [3] . This spline utilized a quadratic distribution 
of the doublet strength over the plane panel obtained by projecting the four corners of the 
curved panel onto the average plane. The six parameters associated with the doublet strength 
allow the doublet strength either to be continuous across the panel edges only in a root mean 
square sense or to be continuous only at panel comers. Discontinuities in doublet strength 
and non-con tiguous panel edges introduce vortex lines which in certain panel configurations 
may introduce large disturbances in the flow. The parabolic arc profile cambered swept wing 
of zero thickness used in figure 39 was analyzed with linearized boundary conditions but 
with different paneling. A portion of the wing paneling near the supersonic leading edge is 
shown in figure 49. The pressure distribution along the 3 strips of panels are plotted in 
figures 50, 51 and 52. The six parameter doublet spline results exhibit large oscillations in 
the pressure coefficient near the tip. These oscillations are seen to be eliminated by using the 
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improved panel method with continuous doublet strength and contiguous panels. 

It should be noted that the poor results from the paneling in figure 49 is over a small 
portion near the leading edge of the same wing shown in figure 39 on which such good results 
were obtained. This indicates a sensitivity of the six parameter spline solution to different 
paneling of the same wing. The limitation of the study to the region in figure 49 was chosen 
because a simple analytic solution for the parabolic cambered wing was easy to obtain. The 
complete cambered swept wing was also computed for the paneling shown schematically in 
figure 53a. AC and CE are special Mach lines along which gradients in the velocity or pressure 
are discontinuous. The pressure in the region ACE was very smooth and in agreement with 
theory. In this region all the panel edges are sonic or supersonic and small discontinuities in 
doublet strength across panel edges do not affect the downwash. The results in ECD were also 
good. In ABC a spike in the pressure occurred at a control point near the point B. This con- 
trol point was found to lie along a Mach cone emanating from an upstream corner point and 
the spike was eliminated when the paneling was readjusted. 

When the paneling in figure 53b was used, unstable results were obtained near the tip of 
the wing. Refining the paneling by an order of magnitude as in figure 49 made the results 
worse as seen in figures 50, 51 and 52. The refined paneling increased the likelyhood that 
disturbances from edge vortices could impinge on control points. The only really effective 
cure was to use contiguous panels with continuous doublet strength. 

To test the sensitivity of the nine parameter doublet spline, the supersonic flow over a 
60° flat plate delta wing at a Mach number of^/lTand angle of attack 5.73° was computed 
using the arbitrary, somewhat pathological, paneling shown in figure 54a. The surface stream- 
line slopes at the panel control points are plotted in figure 54b, and compared with the exact 
linearized solution in reference [23] , page 157. The agreement with theory is good, with most 
of the spread occurring in the panels near the vertex. 

The utilization of combined source and doublet panels on thick wings and fuselages 
eliminates interference of the interior wave reflection with the exterior pressure distribu- 
tion by applying boundary conditions which eliminate the interior perturbation flow. The 
paneling shown in figure 55 was used to compute the flow over the forebody of the B1 bom- 
ber. Use of the 6 coefficient doublet spline yielded smooth results and good agreement with 
experiment when the pressure coefficient was computed by means of the doublet gradient 
as seen in figure 56. The pressure coefficient computed by the aerodynamic influence coeffi- 
cients in figure 57, however, should agree with the results in figure 56, but shows considerable 
departure from the experiment (see ref. [28] ). This deviation is attributed to the influence 
of the line vortices emanating from those panel edges lacking continuity of the doublet 
strength and from noncontiguous edges for which the singularities in the aerodynamic influence 
coefficients from adjacent panels do not cancel each other because of the gap even when the 
doublet strength is continuous. Calculations of pressure coefficient from the aerodynamic 
influence coefficients using the improved contiguous paneling and continuous doublet spline 
agree very well with the experimental results as seen in figure 58. 

The two different panel methods were also applied to the Carlson wing using the paneling 
shown in figure 46. Since the leading and trailing edges are supersonic it is possible to analyze 
the pressures on the wing with sources only on the surface. Since there is no influence of the 
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trailing vortex on the wing because of the supersonic trailing edge, the wake also can be ignored. 
The results for source paneling in figures 59 and 60 indicate an unstable solution resulting 
from the interior Mach wave reflections. From figures 47 and 48 we see that the combined 
source and doublet paneling with the six parameter spline show considerable improvement 
with disturbances from panel edge vortices only occurring in the strip 10 near the tip. The 
use of the continuous 9 parameter spline eliminates all disturbance and yields results in ex- 
cellent agreement with experiment as shown in figures 61 and 62. 

An exact solution of the nonlinear inviscid equations for comparison purposes can be 
generated by using the exact solution for axisymmetric flow past a cone given in references 
[17] to [20] . Any plane which does not pass through the vertex of the cone will intersect 
the shock in an hyperbola (see figures 63 and 64). The shape of the surface was generated by 
computing the streamlines which pass through this hyperbola. The flow properties on this 
surface are the same as if the solid surface itself were immersed in the flow. 

Results obtained by the 9 parameter spline for the surface defined in figure 64 were 
initially very disappointing (see figure 65). The problem is in the approximation of the sharp 
leading edge which has a slope of about 1 degree. When the configuration was taken to be 
the curved surface and the x,y intersecting plane, control points near the leading edges of the 
two networks are displaced in different directions which proved to be more significant for the 
thin edges. 

Much better results were obtained when the configuration was taken to be symmetric 
about the xy intersecting plane (see figure 66). The errors at the top near the trailing edge 
(x = 1, y = 0) are undoubtedly due to inadequate paneling in that vicinity, since the panel 
normals exhibit anomalous variations there. These variations are discussed in a subsequent 
paragraph. Were our interest not confined to the upper surface, symmetry would not be 
available to fix the problem in figure 65. A more general approach is to set up a doublet net- 
work on the mean surface between the plane and the curved surface. This has been done with 
satisfactory results, as shown in figure 67. In this calculation the upper and lower surfaces 
were covered with both sources and doublets, with source strength specified to cancel the 
normal component of the free stream and the potential inside set equal to zero. The normal 
mass flux was set equal to zero on the interior mean surface. 

To study in greater detail the difficulty with thin edges, we considered the flow past a 
sharp-edged delta wing whose lower surface is at zero angle of attack and whose upper sur- 
faces are at a small inclination to the freestream. For the case of a supersonic leading edge, 
the results obtained were poor near the leading edge and got worse downstream (see figure 
68). Analysis of the results revealed a small discontinuity in doublet strength at the leading 
edge junction of the upper and lower surface networks. This was eliminated by replacing the 
= 0 boundary conditions at the control points along the leading edge of one of the net- 
works by W u . n = -U .n. The results improved spectacularly, as shown in figure 69. For 
the case of a subsonic leading edge, the quality of the solution, as evidenced by such data as 
the normal mass flux components on the external surfaces, also was much better when the 
boundary conditions were adjusted to insure continuity of the doublet strength from upper 
to lower surface. However, figure 70 shows that the pressure distribution was just as good 
whether or not the doublet strength was continuous. We attribute this difference in behavior to 
the difference in character between the supersonic and subsonic vortices generated at a doub- 
let strength discontinuity. (These computations were made with the program that did not 
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match doublet strength exactly across network edges. With the later program, these difficulties 
may not have arisen). 

The remedy which worked so well with the delta wings improved the results for the 
conical-flow-field surface, but not to our complete satisfaction. Figure 71 shows the pressure 
distribution at the centers of panels along the leading edge to be quite oscillatory. Much 
better results have been obtained with different paneling of the same surface, and, pending 
further study, we are inclined to ascribe the inaccuracies shown in figure 71 to the high aspect 
ratio of the leading-edge panels. 

Another problem encountered with the concial-flow-field surface was a fall off in pres- 
sure on the triangular panels near the top of the bump. We feel this is due to an anomalous 
variation of the directions of the normals to the panels. However, we now find it is not easy 
to avoid such anomalies in the vicinity of saddle-point-like behavior of the surface being paneled. 
In particular, it is not generally the case that the surface normals converge to the correct value 
at a saddle point as the grid is refined. Unfortunately, results for the pressure distribution are 
quite sensitive to any disparity in surface normal. direction. The alternative panelings shown in 
figure 72 differ considerably in the variation of the x-component of the surface normal over 
the triangular panels, but even the one on the right is sufficiently bad to cause a considerable 
drop off in pressure near the top of the bump (although much less than for the case on the 
left). 


Experience with the delta wing of figures 68 to 70 lead to an improved way of matching 
the doublet strength at network edges. On the downstream side of the common network 
boundary the boundary conditions of the vanishing of the lower peturbation velocity poten- 
tial is replaced by zero mass flux penetration of the upper surface. This dramatically improved 
the delta wing pressures. It also improved the solution to the stream surface of figure 63, 
but the results appear to be adversely affected by the high aspect ratio panels along the 
leading edge. 

12.6 APPLICATION OF THE CONTINUOUS DOUBLET SPLINE WITH 
POTENTIAL BOUNDARY CONDITIONS TO SUPERSONIC FLOW 
OVER WING-BODY CONFIGURATIONS 

The paheling used on the wing body combination described in reference [29] is shown 
in figure 73. The mid wing position with zero angle of attack and yaw was chosen. The flow 
has, therefore, two planes of symmetry and paneling is required for only one quarter of 
the surface. The body consists of an ogive nose terminated by a circular cylinder. The wing has 
a symmetric NACA 65A004 profile. The pressure distribution on the fore-body is shown in 
figure 74. The agreement with the experimental measurements is good. Figure 75 shows plots 
of pressure coefficients at 36%, 53%, and 75% of semispan. The pressure distribution for the 
panel method is smooth and in fairly good agreement with the experimental measurements. 

Over most of the cross section profiles the computed pressures are somewhat higher than the 
experimental values. Near the leading edge the gradient of the pressure coefficient from the 
panel method is steeper than the measured distribution. This difference is due in part to in- 
adequate paneling near the leading edge. The actual airfoil has a blunt leading edge while the 
paneling is necessarily sharp since the leading edge panels must be inclined at angles less than 
the Mach angle. 
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The flow over an arrow wing-body combination using the paneling in figure 76, for which 
no gaps appear between network edges, was computed. The pressure distribution on the wing 
is shown at various span wise cross sections in figures 77 and 78. The calculations are in good 
agreement with the pressure measurements. The introduction of far-field aerodynamic influence 
coefficients reduced the computing time by about 40% over the solutions with all near-field 
aerodynamic influence coefficients. Experimental pressures on the upper surface near the 
leading edge at the 65% and 80% span locations are consistently higher than the theory. This 
is attributed to a leading edge vortex which has not been taken into account in the calculations. 

An analysis was performed on a similar arrow wing-body configuration in figure 79 but 
featuring a twisted wing. The spanwise wing paneling was rearranged to obtain a match with 
the pressure orifice locations on the wind tunnel model. The wing body intersection paneling 
was also refined to match the model closely. The paneling resulted in 983 control points out 
of a possible 1 000. The usual panel boundary conditions of 0 =0 on the interior surface and 
cr = -U . n were applied. The solution required about 600 cpu seconds on the CDC 7600. 

A comparison of wing surface pressure distribution with experiment is shown in figures 
80 and 81 for a Mach number of 1 .70 and 4 degrees angle of attack. Good agreement with 
experiment is shown except at the inboard most station. The leading edge peak in the cal- 
culations is caused by a disturbance at the intersection of the wing at the leading edge and 
the body. The effects of this disturbance are also seen in the calculated body pressures shown 
in figure 82. It is thought that this was caused by poor paneling near the wing-body inter- 
section where one of the subpanels was found to be superinclined. Superinclined panels can- 
not represent solid surfaces. 

12.7 APPLICATION OF SUPERINCLINED PANELS TO CLOSE INLETS FOR 
THE COMPUTATION OF SUPERSONIC FLOW OVER NACELLES 
WITH THE COMBINED SOURCE AND DOUBLET PANEL METHOD 

To test the effectiveness of superinclined panels for closing an inlet, the supersonic flow 
over the nacelle shown in figure 83 was computed using the pilot code. A plane superinclined 
network was used to cover the inlet. The exit area need not be closed as in subsonic flow 
since it lies outside the region of influence for the nacelle. 

Boundary conditions of vanishing perturbation potential and perturbation mass flux 
were applied to the downstream side of the inlet and conventional potential type boundary 
conditions were used on the nacelle surface. To provide a disturbance for the superinclined 
network to cancel, four source panels in the plane 0 = 7r/2 were placed upstream a distance of 
the nacelle radius. To have a solution without superinclined panels for comparison, the source 
distribution on the nacelle was set equal to -U . n and standard mass flux boundary conditions 
were applied to the exterior surface. Figures 84 and 85 show the pressure coefficient along four 
aximuth positions on the nacelle. The two calculations are in good agreement. The flow 
without the superinclined panels appears to have interference from the interior flow due to 
the reflected waves from the lip that is effectively eliminated by the superinclined panels. 

Another test on how effectively the superinclined networks absorb upstream disturbances 
was conducted on the same nacelle. The plane superinclined network was placed inside the 
nacelle at axial station 2.25 (see figure 86). On the upstream portion of the nacelle, zero 
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normal mass flux boundary conditions were prescribed on the upper surface of a source net- 
work. This results in an internal flow to be absorbed by the superinclined network. 

Figure 86 shows the exterior surface pressure distribution on the nacelle as a function of 
axial position for two different paneling of the superinclined network. The 48 panel network 
is seen to be too coarse to absorb adequately the interior reflections. The pressure is oscil- 
latory with a wavelength about equal to the reflected Mach wave pattern. Doubling the panels 
by increasing the number of radial divisions smooths out the pressure distribution. Except for 
the region near the superinclined network, the pressure is in good agreement with LighthiLTs 
theory [25] . 

With the successful development of the superinclined networks, the LES 216 super- 
cruiser configuration, which contains superinclined inlets, was successfully analyzed by the 
advanced panel pilot code. The paneling scheme is shown in figures 87 and 88. A total of 674 
panels divided into 31 networks were used to define the configuration. The model has two 
openings for flow into the engine. The diverter, which splits off the low energy of the boundary 
layer flow from the engine inlet is closed by a superinclined network with the mass flux pre- 
scribed as boundary conditions. The boundary of the inlet opening has edges which are sub- 
inclined to the flow on the upper part and superinclined to the flow on the lower part. The 
upper section of the inlet is covered by a subinclined network and the lower side is covered with 
a superinclined network and the mass flux is specified on both networks. The complete solution 
required 612 cpu seconds on the CDC 7600. Computed results are compared with experimental 
data and with results from a linearized boundary condition method similar to FLEXSTAB 
are shown in figure 89. There is very good agreement between the pilot code calculations 
and the experimental data of reference [31 ] . Note that the results from the linearized code, 
in which the fuselage is represented by a simple body of revolution, are clearly inadequate to 
represent the effects of the canopy and inlet on the wing pressures. Here is an example in 
which the present panel method can more accurately model the aerodynamics of the configura- 
tion than the simpler models currently in use. 

Some waviness is apparent in the pressure distribution on the upper surface for 77 = 0.37. 
This was observed in other pressure distributions and was found to be eliminated by increased 
weighting of the upstream panel center values of the singularity strengths in the least square 
fit. This is explained in Appendix C where a detailed description of the singularity splines is 
presented. 
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13.0 CONCLUSIONS 


Theoretical analysis of the singularities caused by discontinuities in doublet strength 
and in panel geometry, combined with the experience of the flow computations described 
in section 12, indicate that it is important to have panel geometry that contains no gaps and 
a doublet spline system that maintains continuity of doublet strength over the surface. It 
was also found that use of a source panel system alone on closed configurations, such as three 
dimensional wings and bodies, produces interior wave reflections that influence the exterior 
pressure distribution. 

This ultimately lead to the development of the present panel method utilizing combined 
source and doublet panels with potential boundary conditions to eliminate the interior pertur- 
bations. The panel geometry with all contiguous edges was achieved by dividing the basic non- 
planar quadrilateral panel into eight triangular subpanels. Continuity of doublet strength was 
obtained by prescribing a quadratic doublet distribution over each triangular subpanel instead 
of just over the panel projection. The resulting panel method appears to be insensitive to 
choice of paneling provided the paneling is sufficiently fine to describe the solution adequately. 

Superinclined panels, that is, panels which are inclined at angles with the free stream 
greater than the Mach angle, were developed and applied to the closing of inlets on nacelles. 
They were found to be effective in absorbing the incident perturbation flow. With the addi- 
tion of superinclined panels, the present advanced panel method is now capable of solving for 
the flow over complicated wing-body combinations with engine inlets and nacelles. 


Boeing Commercial Airplane Company 
P.O. Box 3707 

Seattle, Washington 98124 
1978 
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APPENDIX A: 

DERIVATION OF BASIC FUNCTIONS AND THEIR DERIVATIVES 


Al. BASIC VARIABLES FOR SUBINCLINED PANELS (SECTIONS 5 TO 8) 
For convenience, we have introduced the following parameters: 

m = (y2 — y 0 / ( x 2 “ x l) X=l/m 

The quantity m is the slope of the panel edge defined by 

x m = ( x “ x l) - (y-yi)/m = 0 


We let 


s m =x m + s/m 
r =y/s^ + z 2 

s z = sx m “ z2 / m 
z m =x m 2 + z 2/m2 


A2. BASIC INTEGRALS FOR SUBINCLINED PANELS 

For the basic integrals we have defined 

I n (r, «-/ t n dt/-v/t(t+2r) 

w n = / s n dsA/s m 2 - r2 = s n ds/R 

Qn = y* s n /r 2 R 

R n = y*s n Rds 

P n = z ^ s n Rds/r2 

s n -/ s n Io( r > s m- r ) ds 


(Al) 


(A2) 


(A3) 


(A4) 
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We shall use the subscript m to denote combinations of functions of the form 


w 


mn 


I 


SmS 1 




ds = x mw n + w n+l/ m 


and similarly for the other functions Q, R, and S. 

A3. RECURSION RELATIONS FOR THE FUNCTIONS- SUBINCLINED PANELS 
The integral I n (r, t) occurs in the first integration of 0 in the form 


/ s m r 

t n dt/ x /t(t+ 2 r) = I n ( r , s m -r) 


To find a recursion formula, we consider 
Integration yields 

t n v/t(t+2r) = r(2n+l)I n (r,t) + (n+l)I n+1 (r,t) 


or 

I n+ ](r,t) = [t n - v /t(t+2r) - r(2n+l)I n (r,t)]/(n+l) 
Substituting the limits 0 to s m -r yields 


I n+1 (r.s m -r) = [ (s m -r) n v^m 2 - r 2 - (2n+l) r I n (r,s m -r)] /(n+1) 
When the first three I n functions are expressed in terms of Iq, we obtain 


(A5) 


(A6) 


(A7) 
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I 


l! = R-rI 0 

1 2 = 3r2l 0 /2+(s m -4r)R/2 

1 3 = - 5r3l 0 /2 + (2s m 2-9s m r + 22r2) R/6 (A8) 


where R =\/s m ~ - r2. 


Following a similar procedure for the w n we write 


A ( s n /rT~^)= ns n - 1 (s m 2- r 2) + s^ /m-s) 
ds V s VSm / / — 9 

V s rn - r 


Integrating and rearranging terms yields 


w n+l 




(n+l) (l - m2) 


[s n R - (2n+l )x m w n /m - n (x m 2 - Z 2) w n _!] (A9) 


Now 


q - f s n ds _ f s n ^ds _ z 2 f s n _ 

J r2 v / Sm 2 - r 2 J ys m 2 - r 2 J 


s n ^ds 


and this results in 


Qn = w n-2 " z2 Qn-2 


(A10) 


The function R n in equation (A4) is more easily expressed in terms of the other 
functions by moving the radical to the denominator. This yields 


R-n w mmn w n+2 z ^ w n 


(All) 
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Expanding w mm Q and using the recursion for W2, we obtain for Rq 

R 0 = ( x m 2 - z2 )w 0 /2 + x m w 1 /2m 

The functions S n are easily expressed in terms of the other functions by integration 
by parts. Hence, using equation (A14), we get 


s n = J s n Io (r. s m -r)ds = J s n log [(s m +R)/(s m -R)] ds 
= [s n+1 Io(r, s m -r) + Qmn+2 ~ w n+l/ m ]/ ( n+1 ) 


Using equation (A 10), we may simplify this expression to 

s n = [s n+1 I 0 + x m w n- z2 Qmn]/(n+l) (A12) 

The P n functions can be expressed in terms of the functions Q n and w n since 

P n = z j* s n Rds/r^ = z J* (s m ^ - r^)ds/Rr 2 

then 


P n “ z (Qmmn “ w n) z ( x mQmn + Qmn+l/ m w n) 


As only the functions Pq and occur, we write these down immediately using the recursion 
formula for the Q n in equation (A 10) 

p 0 = x m (zQmo) + zQm 1 / m " z w 0 


and 

P 1 =x m( z Qml) + z [ x m w l -z 2 Q m0 + w i/m]-zwi 


(A13) 
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A4. ZERO ORDER FUNCTIONS 


Evaluating the integral for Iq yields 


io(> 


r>Sm" r ) = / 


s m" r 


dt 1 . 

y/ t(t+2r) 2 


a m 


a m 


3 


(A 14) 


In this and all subsequent integrals, a form is chosen so that the functions vanish when the 
radical R =\/s m 2 - r 2 vanishes and, hence, the real part is zero for R imaginary. Similarly 


w 0 ( s > s m > r > m ) f 


ds 


v^s m 2 -r 2 


a m 

m 


2y/l - m2 


log 


( s m “ ms ) + R\/l - 
(s m - ms) - R<y/l - m 2 


(A15) 


for Im I < 1 and 


m 


w 0 = 


Vm- 


WQ = 


vT 


...2 


tan' 


tan" 


-1 


ms - s m 


, R >/rn 2 - l y 
s “■ Xst- 


a m 


R yi -a 2 , 


(A 16) 


for Im I > 1. Care must be taken in evaluating all inverse tangents in the aerodynamic influence 
coefficients to select the correct quadrant. The correct signs of the numerator and denominator 
of the argument are essential. In the FORTRAN code the function ATAN2(A, B) is used 
instead of ATAN(A, B). 

It is considerably more complicated to derive Qq and Q\. For this we write 

-±- i-v'T 

s + iz r 2 v 


and obtain 


Iq - Ql - izQo -f- 1 - ^2 _ r 
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Introducing a new variable r = s + iz, the integral Iq becomes 


! Q 


/ dr _ _j 

r-s/a + br + cr^ 2 v / a - 


log 


a + br/2 + R\/a~ 
a + br/2 - R>/a~ 


where a = (x m - iz/m) 2 , b = 2 (x m - iz/m)/m + 2iz, and c - (l - m 2 ) /m 2 . 
Substituting a, b, and c and simplifying yield 

1 


Iq = - 


2 ( x m ~ iz / m ) 


log 


x m ( s m + R) - z2 + i z (vm - R/™) 

x m ( s m - R) - z 2 + iz (y m + R / m ) 


where y m = s - s m /m. Multiplying by x m - iz/m and using the recursion formula for Q n 
leads to 


Qml - w 0 /m - izQ m0 



(A17) 


It is convenient to eliminate all parts of Q m \ which depend only upon endpoints, 
i.e., only upon s and s m , since at endpoints s and s m reduce to y - yj, and x - Xj, i = 1,2. 
We shall show that the function 


j log 


^cp 


L cm 


vanishes in the integration around a contour. Now 

s m sds 


/• s m sas r ds 

Qml- wo/m »/ r2 ^ m2 — -y^f= 


Since s m - s/m = 0 represents the panel edge, then ds/m = ds m along the panel edge; we have 


Qm 1 " w C /m ^ p“) ds ~ 
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In order for the integral to depend only upon its end points we must have 


9 / s m s \ 8/J_\ 

^ s m \r 2 R/ 9s \R/ 

That this is so is easily shown by carrying out the differentiation. Thus, when Q ml is 
multiplied by m or x m the function Q m j may be replaced by 

Qml = w 0 /m 


A similar approach may also be used to simplify zQ m Q. Consider 



zx jYjds 

(s m 2 - s 2 ) R 


= tan X U 

x m K 
= sign (z) tan~* 


"( s ~^ s m) 


for Im I > 1 


mx m R 


I z I (ms - s m )J 


for I ml < 1 


Then 


z Q m o - Qi = zj 


s m ds 

r 2 R 



Substituting x m - s m - s/m and writing ds/m = ds m yield 



(A18) 


(A19) 
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It is easily shown that the line integral on the right hand side vanishes when integrated 
around a closed contour since 


a 

J? 

£ 

= _0_ 

s 

3s m 

_r 2 R (s m 2 -s 2 ) R J 

3s 

_(s m 2 -s 2 )R 


Hence, we may replace zQ m Q by the simpler relation Qj when zQ m Q is not multiplied by 
m or x m , i.e. 


z QmO = Qi 


(A20) 


In evaluating the functions Pq and defined in equations (A13), we are not free to 
use the simpler versions of zQ m Q and Q m \ since these quantities are multiplied by x m and 
1/m and, hence, do not cancel in integration around a panel corner. The complete function 
z Q m o is found by taking the imaginary part of equation (A1 7). For subsonic panel edges, 
m < 1 , we write zQ m Q in the form 


z QmO = " 


1 

2 


tarr 1 


Z (y m - R /m) ~ 
z 2 - x m (s m + R) 


- tan - ' 


z (ym + R / m ) 
z 2 - ( s m " R ) 


For supersonic panel edges it is convenient to combine the inverse tangents. Thus 


zQmO “2 tan 1 


2zR (x m s-z 2 /m) 


2z 2 R 2 - (x m 2 + z 2 /m - z 2 ) (s 2 + z 2 ) 


For points outside the Mach cone, R is set to 0 but inside the Mach wedge of figure 90 
defined by the envelope of Mach cones from points on the supersonic edge for which 

x m 2 + z 2 /m 2 - z 2 >0 

the function zQ m o takes on the values ± tt/2. 

The function Q m j is found by taking the real part of equation (A17). The result can 
be simplified since a common factor 


x m 2 + z 2 / m2 ~ z 2 

can be taken out of the numerator and denominator of the logarithm argument.. The result 
is 

Qmi = wo/m-I 0 (s m , R) 
where Iq is defined in equation (A14). 
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A5. EXPANSIONS OF w 0 AND w x FOR MACH LINE PANEL EDGES 

When Im I = 1 , the panel edge is aligned with a Mach line. The coefficient of the log term 
and arc tangent term in wq then becomes infinite and an expansion is needed for both wq and 
wj. From equation (A 15) we write 

w 0 = (rn/2/J m ) log[(l+0 m z R )/ (l-/3 m z R )] 

where Pm - y/T-ln^ and z R = R/ (s m - ms) . Since j3 m is small, the logarithm term is 
expanded in the form 

wq = mz R 2 £l + z r 2 /3 + z R 4 (l-m2) 2/5 + ZR 6 (i_m2)3/7 + . . 

This expansion is valid for Im I greater or less than 1 . For m = ±1 , it reduces to 

w 0 = R/( s m ± s) 

To obtain a similar expansion for wj we consider equation (A9) and obtain 
WJ =m 2 [R-x m wo/m] / (l-m 2 ) 

Substituting the expansion for wq yields 

wi = m -J •{ R-x m z R [l +z R 2 (l -m 2 )/3 + z R 4 (l -m 2 )2/5 + ZR 6 (1 _ m 2)3/7 + . .] 

1 - ' 

Now 


R “ x m z R ^ 


s m - ms J 

(l -m 2 ) Rs 

m (s m - ms) 


' _ s m - S / m 

s m “ ms 

(l-m 2 )z R s 

in 


then w\ becomes 


W 1 = mz R {s - mx m zR 2 [1/3 + z R 2 (l -m 2 )/5 

+ z R 4 (l - m2)2/7 + Zr 6 _ m 2)3/9 + . . .][ 
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Similar expansions may be derived for w n for n > 1, but these functions are not required 
in the aerodynamic influence coefficients in the form now being used. 

A6. EFFICIENT COMPUTATION OF w 0 , wj AND Qj 

In the calculation of the influence coefficients, the functions w q, wj , and Qj are evaluated 
for each edge at the two endpoints. This involves the utilization of the computer subroutine 
for inverse tangent and logarithm twice for each edge. Considerable speed of computation can 
be achieved by combining the arctangents and logarithms for each edge. For the supersonic 
edge we have 

wo= 7^ ,an '' fcfe) 

i 

where y m = Xs m - s. Using the formula for the tangent of the difference of two angles yields 


1 


WQ = 



tan 1 


/ yi ~^ 2 (y m i R2-y m 2 RQ 
V y m iym2 + C 1 -* 2 ) r 1 r 2 


where the subscripts denote the quantities at the designated endpoints of the edge. Similarly, 
for Qj we have 


Qj = tan ^ 



2 

1 


/ zx m(yml R 2 ~ym 2 R l) \ 
\ x m 2R l R 2 + z2 ym2y m l / 
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For subsonic edges 


w 0 = 


m 


2 (l - m^) 


log 


s m - ms + R-v/1 - m^ 
s m - ms - R \/l - m 2 


2 

1 


m 

(l - m 2 ) 


log (s m - ms + Ry/l- m 2 ) - log (m 2 z m ) 


2 

1 


Since z m (equation (A3)) has the same value at all points of the line, we obtain 


m 

wq = . log 

V 1 - m 2 


y m 2 + V^ ~ m2 r 2 
y m l +v/l -m 2 Rj 


where y m = s m - ms. For Qj we have 


Q\ = sign (z) tan"* 


r£m R l 

2 

U z| ym_ 

1 


= tan - * 


z x m (ymlR2-ym2 R l) ~ 
z2 ymiym2 + £ m 2 RlR2_ 


For edges aligned with the Mach lines, wq becomes singular and we must expand the 
function in powers of _ \2 Letting 


Zr = ymiR2-ym2R] 

ymiym2+(l -X 2 )RjR2 


and expanding the inverse tangent in wq yields 


w 0 = z r 



( 1 - X 2 ) z r 2 
3 


+ 


(l -X 2 ) 2 z r 4 

5 


(l-X 2 )3 Zr 6 

7 + 
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A similar expansion for wi yields 


wi = ~~2 [ Xx m w 0- R 2 + R l ] 


“ j ^2 [ Xx m z r “ r 2 + R l] z r 
Substituting for z r and simplifying yields 



(l - X 2 ) Z r 4 (l - X 2 ) 2 z r 6 
5 + 7 


W 1 = [ ( y m 1 s 2- R l 2 ) R 2-(yml S 1 ~ R 2 2 ) R ]]/[ym2 Vml + 0 “ x2 ) R l R 2] 
_ Zf |V _ (l - X 2 ) z r 4 + (1 -X 2 ) 2 z r 6 _ 1 


The functions wq and wj are seen to exist for the value I X I = I m I = 1 . 

A7. GRADIENTS OF THE BASIC FUNCTIONS FOR SUBINCLINED PANELS 

The functions used in the influence coefficients are functions of s, Sj^, and z when x m = 
s m “ s / m is substituted. When the limits y - y\ and y - y 2 are substituted for s into the 
functions for evaluation at the end points of the panel edge, the quantities s m and s take on 
the values x - x\,y -y\ and x - X 2 , y - y2, respectively. Thus, to find the derivatives with 
respect to x, y and z, we need only compute the derivatives with respect to s, s m , and z. 

Using the operator 


V = (3/3s m , 3/3s, 3/3z) 


we obtain 


[- ms z , m (x m s m z^) , z (s m ms)] 

vwn = r 

m ( z m “ z 2 )R 

Using the recursion relation (A9) to find W] leads to 

m 2 r Rl vx m VWQ 1 

VW1 = VZ^2l~R ~ W °~ ~ J 

where R] = - z) and vx m = (l,-l/m,0) . 


(A21) 


(A22) 
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Using the first of equations (A19), we find that 


VQ 1 = ( Sm 2_ s 2) m ( Zm _ z 2) {<? " s ™ /m ) [ R (- z / x m> z / mx m> 1) 

- z ( s m^~ s > zj) J - (i/ m , -1, o)zr| 


(A23) 


This formula holds for both 0 < I m I < 1 and I m I > 1. 

The gradients of R n and S n can be found in terms of the other functions by differentiating 
under the integral sign, and in the case of y, differentiating upper and lower limits as well. 

This yields 


v R n = ( w mn> R “ w mn /ni, - zw n ) 

v S n = (w n » s n Io(r,s m “ r ) “ w n/ m > - z Qmn) (A24) 


When m = 0, the panel edge is parallel to the panel x axis. If m = the panel edge is 
parallel to the local y axis. In order to represent all panel edge angles, for I m I < 1 , for which 
the panel edge is subsonic, we introduce new variables 

x m = mx m = ms m - s H = ms z = sx m - z 2 

z m = m 2 z m = x m 2 + z 2 (A25) 

and define new w n function by 

w n = w n /m (A26) 


For the supersonic panel edge I m I > 1 and we replace 1/m by X but retain the original form 
of the variables instead of the relations in equations (A25) and (A26). In this form the 
functions hold for all values of m, 0 < I m I < 1 and for all values of X, 0 < I X I < 1 . 

A8. ANALYSIS OF SINGULARITIES OF THE FUNCTIONS FOR SUBINCLINED PANEL 

AND THEIR DERIVATIVES 

A8. 1 SINGULARITIES ON THE MACH CONE 

To find the effect on the flow pattern of discontinuities in doublet distribution and 
source distribution, it is essential to understand the singularities of the basic functions Wj 
and Qj and their derivatives. From equation (A1 5) we see that for subsonic panel edges 

WQ ~ R 
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near R = 0, the Mach cone through an endpoint of a panel edge. The gradient Wq then becomes 
infinite like 1/R as the field point approaches the Mach cone (see fig. 91). This is true for both 
subsonic and supersonic edges. When I m I > 1 , the panel edge lies outside the Mach cone through 
the comer point, wq does not vanish on the Mach cone, however, but letting R -> 0 in equa- 
tions (A 16) yields 


WQ = ±7T/2y/\ - (A27) 

Since x m is always positive downstream of the panel edge, the sign is according to the sign of 
s - Xs m . Using the recursion relation in equation (A9), we see that for subsonic panel edges 

W] R 

near the Mach cone. For the supersonic panel edge 


W 1 = Xx m wo/(x2 - l) = ±Tr\x m /2y/ (l -X2)3 (A28) 


on the Mach cone. The quantities wq and wj are defined inside the region between the 
Mach cones from the two ends of the panel edge and Mach wedge formed by the planes 
containing the edge and tangent to the two end Mach cones (see fig. 92). This wedge region 
is defined by the relation 

z m - > 0 

In this region the value of wq and wj take the limiting value on the Mach cone given by 
equations (A27) and (A28). 

On z = 0, wq when evaluated for each end point is a function which is conical about the 
end point as illustrated in figure 91 . The combination of wq evaluated at the endpoints 

w 0 1 2 = ( w o)l - ( w o)2 

leads to the total contribution from the panel edge described in figure 92. 

The other function which is also defined in the Mach wedge is Qj. The variable x m is 
zero on the vertex of the Mach wedge and is positive within it. Thus, it is seen that for 
I ml> 1, 

Ql = ±tt/2 (A29) 


as R 0 with the sign determined by z (s - Xs m ) while for subsonic panel edges, Qj goes 
to zero like R 0. Since Qj ~ R, then gradient ot Qj becomes infinite like 1/R near the 
Mach cone. 
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The value of Qj for z — 0 in the Mach wedge is easily shown for I m I > 1 . We must con- 
sider the contribution from the endpoints of the panel edge: 

Qn - Q12 

Downstream of the panel edge x m > 0. The quantity s - s m /m is proportional to the oblique 
coordinate to x m , say y m . 

For points x } y lying between 

y m l = si -Smi/rn = 0 


and 

y m 2 = s 2 - s m2 /m = 0 

the quantity y m 2 >0 and y m j <0. For this case Q\\= - 0/2) sign (z) and Qj 2 = 0/2) sign (z) 
and we have 

Qll “ Ql 2 = - sign (z) (A30) 

For points outside for which y m i and y m 2 have the same sign and R = 0, 

Qll - Q12 =0 


For subsonic leading edges and z 0, the function Qj = ±ir/2 if R ¥= 0. The sign is 
determined by the sign of z and of mx m . Consider a point downstream of the intersection 
of two subsonic panel edges as illustrated in figure 93. The quantity mx m is positive for the 
right hand edge in figure 93 and negative for the left hand edge. Since the influence 
coefficient from a panel is found by going around the panel in the counterclockwise direction 
the value of Qj from equation (A19) for the two sides evaluated at the common endpoint 
for I m I < 1 is 


- 7 r sign (z). 


For points between the Mach lines and the nearby panel edge the contribution from the 
function Qj is zero. In the region between the Mach cones inside the Mach wedge for a super- 
sonic panel edge equation (A30) holds. 
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A8.2 SINGULARITIES ON PANEL EDGES 


For z = 0, we have, for subsonic panel edges 


wq 


m 

2y/\ -m 2 


log 


( s m ~ ms) + \/ (l - m 2 )(s m 2 - s 2 ) 
(s m - ms) - >/(l - m 2 ) (s m 2 - s 2 )_ 


Upon multiplying both the numerator and denominator by the denominator, the relation 
takes the form 


w 0 


m 



log 


m (s m - s/m) 

(s m - ms) - y / (l - m2) (s m 2 - s 2) 


(A3 1) 


Since s m - s/m = x m = 0 is the panel edge, we see that wq has a logarithmic singularity at the 
panel edge. From equation (A21) setting z = 0 yields 


vwq 


(-s, s m ,0) 

x m \/ s 2 m - s2 


(A3 2) 


The gradient of wq has a higher order singularity l/x m at the subsonic panel edge. Since this 
formula holds for supersonic edges as well, the gradient has a l/\/s m ± s singularity on the 
two Mach lines through the corner point and is zero between the Mach lines and the panel 
edge. The gradient Qj on the panel is given by setting z = 0 in equation (A23), or 



(s - Sm/m) 

(s m 2 - s 2 ) 


R(0,0,I)j 


dQj/dz = 


(s - s m /m) 
XmVsm 2 - s 2 


9Ql 8Qi 

3x 3y 


(A3 3) 


For I m I < 1 , this has a l/x m singularity near the panel edge x m = 0. For the sup ersonic panel 
edge tm I > 1 , and dQj/0z is zero ahead of the Mach line. Also 3Qj/3z has a l/v/s m ±~s 
singularity at Mach lines. When the panel edge is a Mach line, then the singularity is 1 /\/ x m 
at the panel edge. 
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A9. DERIVATION OF THE BASIC FUNCTIONS FOR SUPERINCLINED PANELS 

(SECTION 9) 

We define the following functions: 


5l ■ xy '/^ and *° = / f 


where 




r 2 = y e 2 + £-z e ) 

2 and R s = >/x 2 - r 2 . 

Let 




t = r 2 = y e 2 + (?-z e ) 2 or f- z e = 

v/t - y e 2 . Then d? = dt/2 y/t - y e 2 

and 


dt 


Or=^ / 


Ql 2 7 

t\/(x 2 -t) (t-y 2 ) 


Integration yields, after some manipulation. 


^ 1 i 

Ql = — tan 1 


x 2 (f-z e ) 2 " ye 2 R s 2 

j_ 2xy e (f-z e ) R s 


Using the formula, tan - * (1/X) = 7r/2 - tan -1 (X), yields 


1 


2Ye R s/ x (^ z e) 

y 2 R s 2 


l- 


2 tan_1 


Qj =-*2 tan~ l 

L J X 2 (f-Z e ) 2 J 
This can be further simplified by applying the formulae 

2 tanfl 


tan 26 = 


I - tan 2 0 


-2y e R/x (f-z)' 


1 " 


y2R2 


L x 2 (f-z e ) 2 J 


and 


tan"l (l/x) = 7r/2 - tan 1 (x) 


Qj = tan * 


x Qe ~ Q ' 

L -y e R s . 


Thus 



The function wq is seen to be 


/x 2 -y e 2 -(f-z e ) 2 


= tan-1 i^e 


The higher order functions w n = f , can be evaluated by a recursion relation. 
Consider J s 


dr - 


nf n_1 Rg -f n (f' z e) 


Integrating yields 


W n +1 = rrr [(2n+l)z e w n + n (x 2 - y e 2 - z e 2 )w n _j - r n R s ] 


from which 


WJ = z e w 0 - R s 


^2 =' ^[Szewi +(x 2 -y e 2 -z e 2 )wo - ?R S ] 
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APPENDIX B: 

DERIVATION OF THE INFLUENCE COEFFICIENTS FOR SUBINCLINED 
PANELS OF SMALL RELATIVE CURVATURE 


Bl. THE CURVED DOUBLET PANEL 

Bl.l EVALUATION OF THE FIRST INTEGRAL IN EQUATION (5.10) TO THE FIRST 
ORDER IN ? 

The first integral of equation (5.10) is expanded in the e’s and only first order terms are 
retained. Using the expansions 


M(x-r + ej -t,yo) = M(x-r-t,yo) + ej n x (x-v- t,yo) 

i = l f 1+ 2e 3- e 2t + ] 

\A (a] t+bj ) % /t(t+2r) L 2 ( t + 2r ) j 



F (x 0 ) dx 0 = 



F( x O)dxo + ei F (x) 


and substituting the new variables 

x m = X - x i -(y-yi)/m, s = y-yo 
s m = x — x j +(yo-yi)/m = X m + s/m 

r 2 = s 2 + z 2 (Bl) 

we obtain 


± 

27 r 



M(x-r-t,y-s) + (z/r) (x-r,y-s) m x (y-r-t,y-s) dtds 

>/t(t+2r) 

M(x-s m ,y~s) ?(x-r,y-s)ds 


z 


2 



M(x-r-t.y-s) [^ xx t-2(^ x + C/r)] dtds 

v/t(t+2r)3 


(B2) 
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To simplify the last integral of equation (B2), we replace the first term in square brackets 
with 


? xx t/(t+2r) = [1 -2r/(t+2r)] £ xx 
and expand the remaining terms to obtain 

T x +?/ r = f x (x-r, y-s) + f (x-r, y-s)/r = fo/ r - tf xx /2 


where f 0 = f( x >y-s)- 
Since 


r^- 

J v / t(t+2r)3 



d(t+2r) 

(t+2r) yj\ (t+2r) [(t+2r) - 2r] 


t 

r \/ t(t+2r) 


integrating the term 


\ C^x + £/ r ) 

v/t(t+2r)3 


of equation (B2) by parts and combining all similar terms in the resulting equation yields 


, a ( r y l r m ~ T 

~2i r 9z | ( 1_Z W 2 ) / 7 

' .r n 


y_y 2 


ju(x-r-t,y-s)dtds 

Vt(t+2r) 


y-: 


+Z 


+Z 


/ 

y~y 2 


K x s m ,y s )l- s m(^o/ r + ^xx/ 2 ) ?Qx 
\/s 2 -r 2 


y-yj 


s m-r 


J ' 1 J ii x (x-r-t,y-s ) 


y-y2 o 


x/ t(t+2r) 


[- 


y - s) + t (^ + 


^xx' 


dtds 


(B3) 
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To perform the integrations with respect to t, we define the functions 


I n < r > t) = 



t n dt 

v/t(t+2r) 


Expanding 

M(x-r-t, y-s) = m - M x t + M xx t 2 /2 (B4) 

where JI = /i(x- r, y-s) 

and F= ?(x-r, y-s) = f 0 “ r ?0x + r2 f xx/ 2 (B5) 

where ?o = 

equation (B3) then becomes 


-L JL 

27 t dz 


0 z fxx/2) 



[m IO" I lM x + l2M xx /2]ds 


+ 2 



M (x-Sm, y-s) [smCiTp/r 2 + r xx /2) - fo x ]ds 
x/ s m 2 - r2 


f /i x f(x-r, y-s) 

L 7 ! 0 

y-Y2 


- 1 1 (m xx ? ( x-r, y-s)/r - M X (ro/r 2 + f xx /2)) 


^2 M xx (?o/ r2 + ?xx/ 2 ) ds 



(B6) 
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where I n = I n (r, s m -r) . 


The integrals in equation (B6) are further simplified by expanding 

M = M(x-r, y-s) and f(x-r, y-s) 

in r and using the recursion relations in equation (A8), we obtain 


2n 8z 


(l-zfxx/2) / [(MO + r 2 fxx/ 4 ) l 0 ' (“Ox " Vs m 2- r 2jds 


y-y 2 


+ z /* M(x-s m , y-s) [s m (fo/ r2 + fxx/ 2 ) - ?0x] ds /\^? " r2 

y-y 2 


- Z 


/* y-y J 

[fOx^Ox + (^0 + 1-2 ^xx/2)Mxx/2] 


y-y 2 


+ z 


/•y-yi s 

y [^OxMxx + (^0 / r2 + fxx/2) (mQx " 2 ^xx)J V s m 2 " r2 


y_y 2 


(B7) 


Since x/ s m 2 - r 2 / r 2 = s m 2 /T 2 y/s^ - r 2 - l/V^m 2 - r 2 

expanding p (x-s m , y -s) in powers of s m and collecting coefficients of 

M(x, y — s), p x (x, y-s), m xx 

yields 


<t>( z) 


= “ 27 T §z y* {m(x, y-s)i//] (s) + p x (x, y-s)i// 2 (s) + M xx t//3 (s)|ds (B8) 


y_y 2 
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where 


'P i - 0 - z W 2 ) *o + 



s m^xx / 2 ~ ^Ox 
\/ s m 2 " r2 _ 


; j. r i /„»- i\ / o ? z [f 0 " s mfOx + s m 2 fxx/ 2 ] 

*2 - - tfOxfo + (zfxx" 1 ) V s m^ - fZ / 

V s m 2 ” f2 


*3 - [(> - 4 - Tl 10 ' * [>0x ‘(>4) t] 



(^0 s m$*0x + s m^fxx/^)/ ^ 


(B9) 


Before performing the integration and differentiation, we now examine the other 
integrals and cast them in a similar form. 

B 1 .2 EVALUATION OF THE SECOND INTEGRAL OF EQUATION (5.10) 

Performing the differentiation with respect to x of the second integral in equation (5.10) 
and introducing the variables in equation (Bl), yields 


000 = - 


_L 

2ir 



f x ( x ~ s nb y~ s ) M (x-s m , y-s) ds 

vV" 2 


27T 




[?xx^( x-r-t > y- s ) 


+ ? x (x-r-t, y-s)M x (x-r-t, y-s)]dtds/ sj t(t+2r) 


(BIO) 
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To perform the integration with respect to t, we expand the quantity in square brackets 
in powers of t. Integration of the second integral in equation (BIO) with respect to t yields 


27 r 


/ 

y-y 2 


y-y i 


C(fxx P + ?0xMx)lQ 


(2f xx M x + fOx^xx)^! + ^fxx^xx^2] 


Applying the recursion relations in equation (A8) leads to 


_L 

27T 


r ^ 

y-Y2 


M + fOxMx) + r ( 2 $xx 


^x + ^Ox^xx) 


+ 9r 2 £ xx M xx /4]lo (2fxxM x + r 0 xMxx)v^^ 


3 (sm-4r) 


fxx^xx r^^ds 


(Bll) 


Equation (BIO) is placed in the form of equation (B8) by expanding \i in powers of r and 
;u(x-s m , y-s) in powers of s m . This yields finally for equation (BIO), with equation (Bll), 


<fiW = - 


_1_ 

27T 



y-yi 

(s) n (x,y-s) + ^2 (s) M x (x,y-s) + ^3 (s) M xx ] ds 


y-Y2 


(B12) 
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where 


n _?x(x-s m >y-s) . 

1_ V^ +Wo 

V s m 2 " 1-2 


n 3 = Sm ?x ^== ? - (?0x " 3?xxs m / 4 ) \/s m 2- r 2 + 3r2f xx l 0 /4 
2v s m 2 ~ r2 


(B13) 


B1 .3 EVALUATION OF THE THIRD INTEGRAL IN EQUATION (5.10) 

Performing the differentiation with respect to y in the third integral of equation (5.10) 
and introducing the variables of equation (Bl), yield 


0(y) = __L f yn s Jy ( x ~ s m . y~ s ) » ( x - s m> y~ s ) ds 


y_y 2 


r \/ s m 2 -r2 


_1 

27 r 


yi _s_ ^ Sm r r xy M(x-r-t, y-s) + f y (x-r-t, y-s)n x (x-r-t, y-s) ds ~[ 

J T J. L >/t(t+2r) J 


y-y 2 0 


1 f y V[ J. f ^ * ?y( x ~ r -t. y)M(x-r-t, yp)dtds 

~ 27r J r i ./C^3 


Sm r 


y-y 2 


Vt(t+2r)3 


Since 


/ 


dt 


\/t(t+2r)3 r-y/t(t+2r) 
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integration of the last integral by parts leads to 


4>(y) = - 


2ir 



ss m [?y ( x ~ s m» y-s)M (x-s m , y-s)]ds 
f2 \/ s m 2 “ 1-2 


2tt 



(r+t)[f xy /u(x-r--t, y-s) 


+ f y (x-r-t, y-s)M x (x-r-t, y-s)] dtds/r 2 y/i t(t+2r) 


(B14) 


To integrate the last term with respect to t, we expand ju(x-r~t, y-s) and f(x-r-t, y-s) 
in powers of t and apply the recursion relations in equation (A8). After considerable simpli- 
fication, equation (B14) then becomes 



y-y 2 


SS 


?y ( x ~ 


'npy 


>m 


, y-s) fi (x-s m , y-s) ds/r 2 v / Sm 2 -r 2 



y-y 2 


) (2fxy^x + fy^xx) * 0 s /2 + ^xy^xx s \/ s m 2~ r 2 


+ fxyF + ?y ( x s m> y s ) (Fx s m^xx/-)] s V^ s m 2_r2 / r2 1 ds 


Finally, 0ty) takes the form 


<p( y) = - 


27T 



[chi(s) ju(x, y-s) +<h 2 (s) M x (x, y-s) +4>3(s) m xx ] ds 


(B15) 
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where 


3*1 s ( s m?y/ r ^ ?xy) / y/ s m^ 


4*2 fxy s ^O ” s fy ( x ~ s m> y s )/ \/ s m^ ^ 


*3 “ - ?y sI 0/ 2 + ?xy S\/s m 2 -r 2 + ss m r y (x-s m , y-s) /2v/s m 2 -r 2 


(B16) 


Dl:4 INTEGRATION OF THE COMPLETE POTENTIAL 

From the form of equations (B8), (B12), and (B15), we see that the velocity potential 
from the doublet distributions is 

0 = <t>( z) + <pW + <t>( y) 


j. f 

2tt 7 


[co^sMx, y-s) + W2(s)At x ( x = y-s)+ w 3 (s)M xx ]ds (B17) 


where coj = -r— z) + £2j(s) +<t>i(s) i = 1,2,3 


(B18) 


Since /lx is a quadratic distribution of the form 


ju = ag + ajx + a2y + a3xy + a4X^ + 


(B19) 


then 0 may be expressed as 


0 2 a n 0n 

n=0 


To separate 0 into its components, we define 


■v-± / 


CO | (s) skis 


(B20) 
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, r n 

J 


co2(s) sMs 


y-Y2 


(B21) 


i r y ~ y] 

x= " 2 ^ / W 3 (s)ds 

y— y 2 


(B22) 


Then the potentials 0 n become 

0 O = <*0 
01 = xaQ + 70 

0 2 = y a o _ a i 

0 3 = xyng - xoq +Y70-71 


0 4 = x 2 a 0 + (2x-x m ) 7Q “ 71 /m + X 

05 = y^aQ - 2ycq + «2 

(B23) 

The basic integrals needed in evaluating a^, 7 j and x 

Rj =/siy Sm 2- r 2 ds 

(B24) 

s i = J sn 0 (r,s m -r)ds 

(B25) 

Qj - J" s i ds/r 2 y / s m 2 - r 2 

(B26) 

Wj = /* s i ds/ v /s m 2 -r 2 

(B27) 
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The functions Rj, S^, Q v and W- are derived in Appendix A. It is also convenient to use the 
subscript m to define functions like 


r Sirius f 

m i J 


sl (x m + s/m) 
s/sm 2 -!" 2 


= x m wi + w i+1 /m 


(B28) 


Similar notation is applied to the functions S, Q, and R when required. 

Using equations (B9), (B 1 3), and (B16), we obtain for equations (B20), (B21), (B22) 

= {’ dz” K 1 ” z ^x/ 2 ) Sj + z (A . Q m j + f xx w m j/2-B* wj)] 

” fxx ($i " w mi) “ B • Wj - fyQmi+1 + fyyQmi+2 + fxy w i+l|‘ 


(B29) 


7i = “ 27 t { ~dz C“ zB * S i + ( z ^xx~0 Rj - z (A . Wj - B • w m j 
+ ? xx w mmi/2) ] + B • Sj - B . w m j + ? xx w mmi " 2? xx^i 
~ ?xy ^i+i - ?y w i+l + fxy w mi+l + £yy w i-*-2 ^ (B30) 


X = -2^{^[(l- 3 % 2L )^ 2 So + S2)/ 2 -zA. s o-zB. s m0 
+ 2zBR 0 -J (l -% L ) R m0 ] + B . S m0 - 2 B . R 0 - f xx R m o / 2 

+ 2~ ^xx( z2 ^mO + ^m2)“fxy^ml “fy^l + fyy^2 + 2 fxy^l|' (B31) 
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where the terms of the form A • Wj and B . wj are defined by 


A . Wj = fWj - fyW i+ i + f yy w i+ 2/2 

(B32) 

B . Wj = ? x Wj - ? X y w i+l 

(B33) 


Since the integrands of Sj and Rj are continuous, differentiating under the integral sign is per- 
missible and we obtain 


9Sj 9Rj 

~dz = - z ^ m i’ ~dz~ ~ ~ ZW i 


(B34) 


With these relations, performing the differentiation with respect to z of the terms indicated 
and simplifying leads to 


a i = 2 tt | z Qmi " ^xx (Sj+z 2 Q m j)/2 

“ ^ * Qmi z ) ~ z (?xx w mi/2 - B i w i) 

+ ?xx w mi/ 2 - fyQmi+1 


+ fyy Qmi +2 + fxy w i+l 


} 


(B35) 


7i = 2?r {" zw i + ^xx ( z2w i + R i) + B * ( w mi ' z2 Qmi) 

+ (l + z jj“) [A • Wj - B • w mi + ?xx w mmi/ 2 ] + ^xx w mmi 

+ ?xy s i+l + fy w i+l " ?xy w mi+l " fyy w i+2 | (B36) 


X = - 2^ { zS 0 + (fxx/ 2 )[ z2w mO " 3z 2 S 0 /2 + 3S 2 /2 - R m o/ 2 ] + z2b * (QmmO " 2w o) 
- fxy^ml “ ?y B l + ?yy^2 + 2 ^xy B l + ^ • (z 2 Q m o - Sq) ^ 


(B37) 
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When curvature is neglected, f = A . = B . = 0 and 


a i = z Qmi/ 2ir 

7i = - ZWj/27T 

X = - zS 0 /2tt (B38) 


the potential then reduces to equation (5.17). 

By examining the terms in the velocity potential involving f , we find that the potential 
has the same singularities as the velocity components for the flat panel. From equation (B35) 
for cq), we find the term 

(A . zQ m0 ) = -^-(fQl - zQ m i ? y + Q m2 ? yy /2) (B39) 

leading to the derivative 


3Ql/3z 

From equation (A23) we see that the gradient of Qj has a 1/R singularity on the Mach cone 
R = 0. The velocity components then will have 1/R 3 singularities on the Mach cone. This 
singularity is too strong since it is not integrable and renders the terms involving f impractical. 
Expanding the potential to the first order in f acts like a differentiation of the flat plate 
potential. To eliminate this singularity, the doublet distribution would have to be prescribed 
on the curved surface itself. It was decided that flat panels only would be used. 

B2. THE CURVED SOURCE PANEL 

We shall evaluate the influence coefficients for the source panel in a manner similar to 
that followed for the first integral of the doublet panel in the preceding section of B 1.1. 
Equation (8.2) differs from the first integral of equation (5.10) only in the differentiation 
with respect to z. Consequently, expanding this integral, retaining only first order terms in 
J can be written down from equation (B3). This is 


0 = - 


_L 

27T 


{(1 - zfxx/2) f 


y-yi 


y_y 2 



a(x-r-t,y-s)dtds 
y/ t(t+2r) 


+ z 


/ 

y-y 2 


y-yi 


p ( x s rri’y s ) H s m(?o/ r2+ £xx/ 2 ). ~£qx] ds 
V /s m 2 - 1 ’ 2 
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+ z 


y yi r Sm r g x & (x-r,y-s) It + t(?p/r 2 +? xx /2)]dtds 


If 

y-y2 o 


y/i t(t+2r) 


where fg = ?(x> y_ s) and s m and r are defined in equations (A3). Now 

a(x-r-t,y-s) = a(x-r,y-s) - a x t 

since a is linear in x and y. Performing the integration with respect to t yields 

/* y - y l , 

<p = - 57 J Jo- Z W 2 ) [ff(x-r.y-s) Ig(r,s m -r) 
y-y 2 

o (x-s m ,y-s) [s m (f 0 /r 2 + $ xx /2) - f 0x ] 

IlO=Sm-r)J + z y== 

v s m Z “ r 


-o. 


+ za* 


ifatL) lQ (r , Sm _ r) + (2f + I , (r,s m -r) 


ds 


Eliminating Ij by means of equation (A8) and expanding 

cr(x-r,y-s) = a(x,y-s) - ro x = oq - r o x 

lead to 

y-yi 


0 


27J- f { 0 ~ z fxx/ 2 ) [ao ^0 ^xR] 
y-y 2 

+ za (x-s m ,y-s) [s m (fg/r 2 +f xx /2)-fg x ]/R 


+ za, 


,[-foxio + (^| + % L ) R ]} ds 


where 


R 


-v^T2. 


(B40) 


(B41) 


(B42) 


(B43) 
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To express the integrals in terms of the basic functions of appendix A, we substitute 


R s m 2 1 

r 2 t 2r R 

We then expand a(x,y-s) and o (x-s m ,y-s) in s and s m . After combining and simplifying we 
obtain 

1 f y ~ Yl 

0 = - 2 ” J [a(x,y)0 1 (s) + a x 0 2 (s) + o y 03 ( s )] ds (B44) 

y-yi 


where 

0 1 = (l-zf xx /2) Iq + zs m fo/ r ^ + z ( s m?xx/2 ~ fOx) /R 

*2 --*♦*» 

03 = - s (l-zf xx /2) Iq - z s m £*o s/r^R - zs (s m f xx /2 - £q x )/R (B45) 


where = ?( x >y-s) 

Since the source distribution is linear in x and y 

o = aQ + ajx + a 2 y 

then the velocity potential may be expressed in the form 

0 = a O0O + a 101 + a 202 


Then 


00 = * 


_L 

27 r 


/ 

y_y 2 


y-yi 


^ i (s) ds 


(B46) 
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i r y ~ yi ( * 

~ ~ 2i t J { 0~ z fxx/ 2 ) ^0 + zs mfo/ r ^^ + z (fxx s m/ 2 ” ?0x) /Ryds 


y-y 2 


and integrating yields 


^0" 27[ s o + z xo] 


where 


XO A • QmO + ^xx w m0/ 2 - B • wq - f X x^o/ 2 


and A . and B • are defined in equations (B33) and (B34), Now 


i r~ yi 

^1= x 0O-27 / *2(®) ds 

y— y 2 


Finally integration yields 


where 


4>l = x0Q + Ro/ 2,r + z Xi/2tt 
Xj = A • w 0 - B • w m0 - ?xx w mmO/ 2 


Similarly, 


I 


/ 


y-yi 


0 2 = y*o " 27 / ^3 (s) ds 

y-y2 


(B47) 

(B48) 

(B49) 

(B50) 
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and 


02 = y0O + Sj/2ir + zX2/ 2ir 


(B51) 


where 


*2 = fxx(w m i - Sj) /2 + A • Q m i - B . wj (B52) 

To compute the perturbation velocity vector, it is convenient to write the complete 
velocity potential which becomes 


0 = - (So + zxo) - CT x (Ro + Z X l) - °y (S 1 + ZX2) | 


(B53) 


In the same manner as for the doublet panel method, we obtain the relations for the velocity 
by the gradient operator 


'JL JL A\ 

,8s m ’3s ’ 3z / 


Thus, we obtain 


v 0 = - ^ r {a(x,y)[vS 0 + v(zxo)] 

+ (ff x ,a y , 0) (Sq + zxo/2tt) 

-a x [v R 0 +v(zxi| 

-a y [vSj +v(zx2)]| (B54) 


Since the expression for the velocity potential contains the term 

z QmOf 

in zxo °f equations (B53) and (B48), we see that the contribution to the source velocity 
components from the panel shape parameters have some of the same singularities as the 
flat doublet panel. As we discovered from the doublet, expanding in the shape parameter f 
acts as a differentiation of the basic flat panel aerodynamic inlfuence coefficients. 
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APPENDIX C: 

DESCRIPTION OF THE SPLINES FOR SUPERSONIC SOURCE 
AND DOUBLET PANELS 


Cl. SOURCE PANEL 

For each panel, an average plane is defined which passes through the panel center and 
midpoints of the lines joining the panel corners. On the panel projected onto this average 
plane, the source strength a is taken to be linear in the coordinates of this plane, i.e., 

a = erg + a x x + cry y (Cl) 

The three parameters a Q , a x , and o y which are needed to specify the source strength on each 
panel are not, per se, basic unknowns of the problem. Rather, they are calculated as needed 
in terms of the source strengths at the centers of the panels, which, in turn, are determined 
so that the solution meets the specified boundary conditions. 

The parameters of equation (Cl) are determined by a weighted least-squares fit to the 
source strength at the centers of the eight surrounding panels. That is, if the source strength 
at the center f j, 77 j of the ith of the 9 panels involved is called cq the quantity 

9 

E = 2 w i (°i ~ °0 ~ °x % i - °yWl)^ (C2) 

i = 1 

is made stationary with respect to variations in a Q , a x , and Oy. The weight Wj is equal to 
108 if i is the label of the center of the panel for which equation (Cl) is sought, and Wj = 1 
otherwise. 

If the panel is at the corner of the network, there are three neighboring panels, which 
is enough to make this system determinate (see figure 94). 

Two points are worth noting about this representation of the source strength. First, 
there is little reason for the source strength to be continuous from one panel to the next, 
since source discontinuities do not introduce infinite singularities in the flow as doublet 
strength discontinuities do. Of course, in the limit as the grid is refined, we expect the 
discontinuities in the source strength to vanish, so that these discontinuities furnish a 
measure of whether the mesh is sufficiently fine. Secondly, since the source distribution within 
the panels is determined in terms of the source strength at the centers of nearby panels, the num- 
ber of unknowns governing the source strength is exactly equal to the number of panels. 
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C2. DOUBLET PANEL 


Since it was shown in section 3.7 that the jump in tangential velocity across a singularity 
distribution is related to the local value of the derivatives of the doublet strength, while the 
jump in normal flux is given in terms of the source strength, it is necessary to use a higher 
order representation of the doublet strength than of the source strength to acliieve the same 
level of representation of their velocity fields. Thus, the doublet strength is represented by 
a quadratic: 


ju = Mo + Mx x + M y y + Mxx x2 / 2 + M X y x y + My y y 2 /2 


<C3) 


For the six parameter spline, a quantity similar to equation (C2) is formed for the 
doublet strength in equation (C3), namely, 


E = 


N 

E 

k-1 


w k (^0 + ^x£k + *V?k + ^xx£k 2 / 2 + ^xy£k ^k + Myy% 2 / 2 " \) 


2 


(C4) 


where is the value of the doublet strength at the control point of the kth neighbor- 

ing panel (or point on network edge) projected upon the average plane of the panel for which 
the coefficients are to be defined. The quantity E is made stationary with respect to the 
coefficients of the quadratic defining the doublet strength. To facilitate the determination of 
these coefficients, we define the vectors: 


c i “ (M0» Mx> My 5 Mxx? M\y> Myy) > i 1, 2, 3, . . 6 

a ki = 0 - £k> ’ik’ sk 2 / 2 ’ ^k^k’ ^k 2 / 2 ) k not summed 

1 = 1 , 2,. ..,6 


(C5) 


then 


N 

E= Z ( c i a ki “ ^k) 2w k 

k= 1 


from which 


3E 


z 


k = 


w k ( c i a ki “ ^k) a jk 0 
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and 


N ) N 

k a kj a ki w k j c i = k J j w k^k a jk 


Solving for Cj yields 


( N ] 

- 1 ( N ) 

c i - \ 2 a kj a ki w k i 

\ 2 w k X k a jk v 

U« i J 

\ l k = 1 ) 


(C6) 


Equation (C6) may be expressed as the product of an N by 6 matrix with the vector of 
parameters X^: 

(on ) N X 6 

c i“ j^ikj Xfc ( 

This is the basic relation for evaluating the coefficients for the six parameter spline. 

For the nine parameter spline, a relation similar to equation (C 4) is used but the spline 
is made to fit exactly the value of X^ which corresponds to the center value of the panel for 
which we are seeking the doublet distribution. A weighted least squares fit is then made with 
the values of X for the adjoining panels. It was found that smoother pressure distributions 
were obtained when the weights for the upstream panels were increased, with the highest 
weights (about 7) for the panels directly upstream. This weighting recognizes the predomin- 
ance of upstream influences in supersonic flows. 


As was shown in section 6, it is essential that the doublet strength be continuous across 
all panel edges if we are to avoid the infinite singularities emanating into the flow field from 
panel edges. This is accomplished by dividing the panel into eight subpanels and defining a 
quadratic distribution over each subpanel. The determination is accomplished in the stages 
following. 

For each panel, nine doublet strength parameters are defined as the values of the 
doublet strength at the nine vertices of the eight triangular subpanels strength at the nine 
vertices of the eight triangular subpanels by using equations (C3) and (C7). Since the 
definitions of 8 of these doublet strength parameters are associated with points on the' 
boundaries between adjacent panels, they are exactly the same for the adjoining panels. It 
will be shown that this achieves the desired continuity of doublet strength across panel 
boundaries for a quadratic doublet distribution. 

Six of the ten straight lines which bound subpanels contain three subpanel vertices. In 
figure 96 these are the lines ABC, DEF, GHI, ADG, BEH, and CFI. The doublet strength 
parameters associated with the vertices may, therefore, be fitted with quadratics along each 
of these six lines. 
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Since three doublet strength values along each line uniquely determine a quadratic, the 
doublet strength on a panel edge not only agrees with the doublet strength on the adjoining 
edge at the three common points, but also at all other points of the line. By using piecewise 
quadratic distributions in the eight subpanels of figure 95, we can construct a nine parameter 
spline which is continuous everywhere in a network. This can be easily seen by considering 
figure 97, which shows values of doublet strength interpolated by the quadratic at the mid- 
points of the subpanels edges. Every subpanel has six values of the doublet strength which 
uniquely defines the quadratic inside the triangle. As seen in figure 98, continuity across 
subpanel edges is assured since the resulting quadratics from adjacent subpanels are uniquely 
defined on the panel edges by the three common values of doublet strength. 

In place of actually determining the doublet strength at subpanel midpoints, a quantity 
is defined for each edge which provides the required additional information. This quantity 
is 


Kij=M(Pi)+ VM(Pi).(Pj-Pi)/2 

= J“(Pj) + V/i(Pj) * (^i " Pj)/ 2 


(C8) 


where Pj and Pj are the position vectors of the two endpoints of the subpanel edge. For a 
quadratic distribution, these two forms can be shown to be equal. Using equation (C3) for ju 
with points r)- x and $[, we obtain 


K ij = Mq + M x (£i+£j )/ 2 + M y Oli+i?j) / 2 

+ Mxx(£i£j)/2 + M X y(^i Vj + %i 1 ?i)/2 
Myy (C9) 


which is seen to be symmetric in i and j. These quantities are obtained on the 12 subpanel 
edges corresponding to the lines DEF, GHI, ADG, BEH, and CFI by using equations (C3) 
and with the coefficients defined by equation (C7). Since the gradient of a piecewise quadratic 
is not necessarily continuous everywhere, we approximate the Ky value for the edges DH, HF, 
FB, and BD by a linear combination of the Ky for surrounding edges. For example, the value 
of Ky for the side DH is found by a linear combination of the Ky for the sides DG, GH, 

DE, and HE. For this evaluation we need to develop parameters associated with the departure 
of the panel from a parallelogram. 

We now calculate parameters^ which measure the skewness of the panel. Referring 
to figure 98 we define the vector P mn = P m - P n , where P m and P n are the position vectors 
of the comer points. Then, for the upper quadrilateral of the panel, P1P6P9P5 in figure 98, 
we define constants C\ i such that 

Pl9 = (l+Cil)P59+(l+C 2 1 )P69 


137 



where the tilda denotes the projection of the point on the average plane. If Cj 1 = C 2 * = 0, 
then the quadrilateral is a parallelogram. We eliminate C 2 by performing the cross product of 
the equation with and obtain 

P 69 X(Pi9-P59) =C 1 1(P 69 XP 59 ) 


Now P ]9 = P 59 + P] 5 , and since the crossproducts are vectors normal to the plane (or 
effectively scdars), we obtain for C\ 1 

_ 1 _ ^69 xP 15 
Cl — C* ~ 

p 69 x ?5 9 

Now the tilda marked vectors differ from the complete vectors by components normal to the 
average plane, hence, the value of C\ * can be expressed in terms of the points Pj by 

c 1 _ ( p 69 xPis) ■ n 
1 ( p 69 x p 59) * " 

since the contributions to the numerator and denominator from the normal components of the 
vectors vanish. Similarly, eliminating Ci * by the cross product of the equation with Peg 
yields 


c 2* =( p 16 xP 59)/(P69 xP 59) 


since P ]9 = + Pi 5 . For the remaining comers, we take the points P 2 P 6 P 9 P 7 , P 3 P 8 P 9 P 7 , 

and P 4 P 8 P 9 P 5 and obtain the following formulas for CjJ : 


c 1 = ( p 69 xP 15) -n 
1 ( p 69 x p 59) * n 

c 2 = ( p 69 xP 27) 

1 ( p 69 x p 79) * " 

c 3 = ( p 89 x p 37 ) • n 
1 (Pg 9 x P7 9 ) . n 

c 4 = ( p 89 xP 45) «n 

1 (Pg 9 x P5 9 ) . n 


(CIO) 
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for C^, 


c i = (Pl6 x ^59) ■" 

2 0*69 x ?59) • « 

c 2 = (^26 x ^79) ■» 

2 (P 69 x P79) . n 

c 3 = (P38 XP79) . n 

2 0*89 x ^ 79 ) •« 

4 = (^48X^59) *n 

2 (^89 x ^59) •“ (C11) 


We are now able to define the Ky quantities for the subpanel edges DH, HF, FB, and BD 
in figure 96. Let Fj = 1 + C\* + C2* . For the value of K53 and side P5 - Pg in figure 98, we 
have 

K58 = TT ?7 ( K 59 + K 89 ) + YTp-j (K51 +K 81 ) 


Similar results may be written down for the other edges, P5 - P5, P5 - P7, and P7 - Pg. 

For each triangular subpanel we now have the values of the doublet strength at the 
comer points and values of Ky for each of the sides. The coefficients of the quadratic in 
equation (C3) then can be expressed in terms of these six quantities by 
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or 



(C12) 


where /z j , (jl 2, M3 are the comer values of the doublet strength on the triangular subpanels, 
and K] , K2, K3 are the values of the parameter Ky for the three sides of the subpanel. 

The matrix D is given by 



1 

1 

1 

1 

1 

1 


£1 

V\ 

£l 2 /2 

£2 

7?2 

£ 2 2 / 2 

£3 

773 

£ 3 2 / 2 

£1 +£2 

VI +V 2 

£l £2 

2 

2 

2 

£2 + £3 

V 2 + V 3 

£2 £3 

2 

2 

2 

£j +£3 

VI +V 3 

£3 £l 

2 

2 

2 


£l *71 
h V2 
£3 173 



j/SlTO \ V\V 2 
2 \H2ni) 2 

J_/ %2V3\ 32211 
2 \H?2£3/ 2 


]_/£l V 3 \ VI 33 
2 \^3 r n) 2 


(Cl 3) 


The matrix D may be inverted whenever the triangular subpanel does not degenerate into a 
line, for which we do not need the coefficients for m* 

Hence, within each panel we have a continuous piecewise quadratic representation of 
the doublet strength in terms of the nine doublet strength parameters. The doublet 
strength is also continuous with the doublet strength on neighboring panels. 

We shall now explain how the nine values of the doublet strength at the vertices of 
the eight panel subpanels in figure 96 are obtained. As in the case of the source distribution, 
it is convenient to regard as the basic unknowns the values of the doublet strength at the panel 
centers. This helps to assure a correspondence between the number of unknowns and the 
number of boundary conditions. For panels well within the network interior, the doublet 
strength parameter at a corner of a panel is determined by least-squares fitting of the 
quadratic in equation (C3) to the doublet strength at the centers of 12 nearby panels, as 
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shown in figure 99a by a means of equation (C7). The squared errors at the four panel 
mispoints closest to the corner point in question are given a weight of 10^ while the outer 
uncircled points are weighted according to a formula which gives greater weight to points 
upstream of the point where the doublet strength is desired, namely 

W k = 1 + M (l - R k • 6) 

where R k is the vector distance between the central point and the kth outer point and C is 
the unit vector in the compressible direction (usually the freestream direction). 

Similarly, the doublet strength parameters at the midpoint of a panel edge is found by 
least square fitting a quadratic to the doublet strength at the centers of 12 nearby panels, 
the squared errors at the centers of the two panels bordering the edge being weighted at 
1 03 while the other points are weighted according to the formula. These arrangements are 
illustrated in figure 99a and 99b. 

The doublet strengths at panel corner points along network edges are required for the 
nine parameter spline. Because of the necessity of matching doublet strength across network 
edges, the determination of the doublet strength at panel comers on a network edge must be 
confined to the values of the doublet parameters along the edge. We recall from figure 14 
that singularities are provided at the midpoints of the panel edges bordering networks and at 
network corners. The corner values are the actual doublet strength. The midpoint edge values 
are the doublet strength for the six parameter spline; but for the nine parameter spline we use, 
instead of the doublet strength, the quantity 

X i = M (Pi-l) + \ VM (Pi-l) • (Pi - Pi-l) 

= M(Pi)4MPi)*(Pi-l-Pi) (Cl 4) 

In terms of these parameters, the comer values of doublet strength are given by 



Pj-Mj-il \ i+i + _ | Mj - Pj | Xj 
| Pi - Mj_i | + | Mj - Pj| | Pj - Mj_i| + | Mi -Pi 


(C15) 


where Pj denotes the comerpoint and Mj the midpoints as shown in figure 100. This can be 
shown easily by substituting Aj+j by the first of equations (Cl 4) and Aj by the second, while 
noting that the gradient terms cancel. 


141 


Equation (Cl 5) is equivalent to straightening the edge and using a single variable to 
represent /i along the edge. In terms of a single variable, then equation (Cl 4) takes the form 


M + l-Mi + i+Y^) ( x i - x i+l) 


(C 16) 


Since we require n at the midpoint, we let xj+j and xj be ±1, respectively. Then along the 
panel edge, the doublet strength varies as 


M(x) = M m +- 


Mi+1 - Mi 


' Mj+1 + Mj 


-Mm ^ 


(C17) 


where denotes the value of the doublet strength at the midpoint of the panel edge. From 
equation (Cl 7), we have 



3u i+l + Mj-l 

2 


- 3 Mm 


which gives for equation (Cl 6), after solving for /i m , 

Mm = Mj+l/4 + Mi/4 + Xj+j/2 (Cl 8) 

Thus, the panel comer values and panel midpoint values of the doublet strength are determined 
by the singularity parameters Xj at the panel midpoints on the network edge by equations 
(C15) and (C18). 

With the values of the doublet strength at panel corners of the network edge defined, 
we are now able to compute the comer and midpoint values for all the panels in the network 
for application to the nine parameter spline. For the panels near the network edges, we 
obtain the corner values of the doublet strength by using equation (C3) and (C7) to obtain a 
least squares fit with neighboring panel center values and with network edge values. 

Figure 99 shows the fit for interior panels well within the network. For panels close to the 
network edge, the panel comer and edge midpoint values are computed with a least square fit 
with neighboring doublet strength values which also include corner points lying on the network 
edge. Figures 101a and 10 Id show a fit made with only ten points. However, the doublet 
strengths at the two corner panel points lying on the network edge are expressed in terms of 
the two basic parameters on the network edge at the panel midpoints of the adjoining panels. 
Thus, the vector of coefficients for equation (C3) is computed in terms of a matrix product 
of a 12 by 6 matrix and the neighboring Xj^ basic singularity parameters. In figures 101b and 
101c a least squares fit of eleven points is used. All cases result in a 1 2 by 6 matrix with the 
1 2 basic singularity parameters. 


142 



Naturally, special attention must be given in cases for which the panels are triangular 
rather than quadrilateral, if only because it is impossible to fit a quadratic to three values at 
a common point. However, the same basic selection of unknowns serves to define enough 
doublet strength parameters to fix piecewise continuous quadratics in each of the six subpanels 
of any triangular panel (see fig. 102). 
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APPENDIX D: 

DESCRIPTION OF THE PILOT CODE FOR COMPUTING LINEARIZED 
SUPERSONIC FLOW OVER AIRCRAFT CONFIGURATIONS 


The basic formulas for the computation of supersonic flow by the improved higher order 
panel method have been derived in the foregoing sections. These derivations are detailed 
and somewhat lengthy with the formulas actually to be computed scattered throughout the 
document. To present an overall view of the method we shall describe, in some detail, the 
basic calculations used in the pilot code to compute linearized supersonic flow and refer to 
the equations and the sections in the text which are applicable. This section is not meant to 
be a description of inputs and procedures to be applied by the user in running the pilot 
code, since a document written for this purpose is available in reference [32] . 

Section 3 presents the basic theory of linearized supersonic flow following to some 
degree the work of G.N. Ward [10] . Section 4 gives a brief general description of the panel 
method. The aerodynamic influence coefficients are derived and a discussion of their 
properties is given in sections 5 through 9 for both the subinclined and superinclined panels. 
Some of the pure mathematical details have been relegated to the appendices A and B. Since 
the formulas for the influence coefficients were derived for subinclined panels parallel to the 
freestream direction and superinclined panels normal to the freestream, the influence 
coefficients must be transformed for panels at other inclinations by means of the transforma- 
tion matrices derived in section 10. If panels are located well within the upstream Mach 
cone from a control point, then the computations can be made more economical by using the 
far field expansions described in section 1 1 . 

Dl. BASIC GEOMETRY AND PANEL DEFINING QUANTITIES 

The first step in using a panel method is to divide the configuration into networks. The 
network boundaries are distinguished by discontinuities in the configuration geometry such 
as wing and body junctions. The spline representation of source and doublet assumes smoothly 
varying geometry so the body is divided into networks and appropriate matching conditions 
are used across network boundaries for the doublet strength. Networks also may be dis- 
tinguished by the boundary conditions applied to them. For this reason, no network can con- 
tain both subinclined and superinclined panels. 

Each network is divided into panels by a double array of grid points on the surface. The 
rows and columns are chosen so that the cross product of a vector in the column direction 
with the vector in the row direction is pointed in the direction of the positive outward normal 
to the surface into the fluid. 

The boundary conditions are applied at the center of the panel and for this we require 
the unit normal vector. The panel normal is chosen as the normal to the average plane or to 
the central parallelogram of figure 98. This is defined as the cross product of the vectors 
formed by the lines joining the midpoints of the sides of the parallelogram or 
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^ ( P 6 +P 7 _ P5 +P 8) X (P6 + P5 _ PTIPS) 

| p6 +p 7 _ P S +P 8 ^ X ^P 6 xP 5 _ p 7 +p 8^ j 

Normals are also defined for the subpanels. Although these are not used for boundary 
conditions, they are used to define the transformation matrices required for the aerodynamic 
influence coefficients of each subpanel. Also computed is the center of each panel or 
subpanel. It is defined as the average of the corner points and is the origin of the coordinates 
in which the formulas for the aerodynamic influence coefficients are expressed. For 
supersonic flow, one of these coordinates is aligned with the projection of the freestream 
velocity on the panel. 

The control points at which boundary conditions are applied are then computed, the 
locations of which are illustrated in figures 13 and 14 for the source and doublet networks with 
analysis boundary conditions. For the nine parameter spline, the panel center control point 
must be away from exact center because of this point being a common vertex to the four 
central subpanels of figure 98. 

The control points at the edges of doublet networks are used for matching doublet 
strength across network boundaries. For the 6 parameter spline and for earlier versions of 
the nine-parameter spline, the edge matching control points are moved into the network an 
infinitesimal distance. This is essential since the infinite singularity at the panel edge in section 
6 is used in the six parameter spline to match the doublet strength at network boundaries. 

For the nine parameter spline, the continuity of the doublet strength allows us to drop all 
line vortices at panel edges and in the earlier version, a simulated incompressible vortex is 
introduced at network edges as described in section 4 with equations (4.1) and (4.2) to match 
doublet strength. The latest version matches doublet strength exactly by prescribing regular 
downwash boundary conditions at the edge control points for the upstream network to 
determine fhe doublet strength and using these values in the downstream network, to obtain 
exact matching of the doublet strength across the two networks. 

Up to this point, the calculations are carried out without regard to whether the flow is 
subsonic or supersonic or for any choice of Mach number. In summary, we have computed 
the following quantities: 

1 . Designation of networks and the double array of panel comer points on the surface 
corresponding to each network. 

2. Additional points of each panel used to define the eight triangular subpanels. 

3. Unit normals to each panel and subpanel. 

4. Location of the control points at panel centers for the boundary conditions and control 
points at network edges for matching doublet strength. 
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D2. TRANSFORMATION MATRICES 


The aerodynamic influence coefficients were developed for a scaled coordinate system 
with subinclined panels parallel to the free stream and superinclined panels normal. To calcu- 
late the transformations which convert the aerodynamic influence coefficients for any general 
inclined panel, we must introduce the Mach number at which the pressures and velocities 
are desired. 

These transformations were derived in section 10. With Xp, yp, Zp designating the 
coordinate system on the subpanel whose origin is at the subpanel center, the transformation 
from a reference system x r , y r , z r to the panel system is described by the equations (10.37) 
(10.38), and (10.39) or 


(V y P’ z p) = i A 2^ M I A rl 


x r " x o\ 
y r - y o 
z r ” z o/ 


(Dl) 


where Xq, yg, zq are the coordinates of the subpanel center in reference coordinates. 

For source influence coefficients there is a correction required which is given by the 
Jacobian of the transformation from the reference system to the panel coordinate system in 
equations (10.42) for subinclined panels and (10.44) for superinclined panels. 

D3. CALCULATION OF AERODYNAMIC INFLUENCE COEFFICIENTS 

The influence coefficients to be calculated depend upon the boundary conditions 
that are to be applied to the surface. The pilot code is set up on the assumption of a 
combined source and doublet panel system. A single singularity surface is selected by set- 
ting the unrequired singularity equal to zero. The most likely boundary conditions are the 
potential type boundary conditions, in which the source strength o is defined to cancel the 
the normal component of the free stream velocity on the outer surface, i.e., 

a = -U« n (D2) 

and the perturbation potential on the interior surface of the configuration is set equal to 
zero, or 


tfr= 0 


(D3) 


These boundary conditions were shown in section 3.8 to satisfy the vanishing.of the normal 
component of the linear mass flux on the outer surface and to eliminate the interior flow 
perturbations. Equation (D3) leads to a set of simultaneous linear equations to be solved 
for the doublet parameters. 
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Since the source parameters are defined by equation (D2) they contribute to the right 
hand side of the simultaneous equations. 

The boundary conditions in equation (D3) is applied to the downstream network edge. 

On the upstream network edge of the adjoining network we apply the boundary condition 
in equation (4.1). 

The contribution from a panel edge to the potential from a quadratic distribution of 
the doublet strength is given in equation (5.17) where the functions are derived in Appendix 
A. The contribution of the entire subpanel to the perturbation potential at the control point 
is found by evaluating and summing equation (5.17) for each end point of the side and for. 
each side by proceding counter-clockwise around the panel. The quantity jtx and its derivatives 
are defined in equation (C3) evaluated for the x, y, z coordinates of the control point where 
the boundary conditions are to be computed. Substituting equation (C3) for fi into equation 
(5.17) and combining coefficients yields a relation of the form 

0 = Aj c j (D4) 

where the vector Cj is defined in equation (C5). The contribution to the velocity potential 
from the panel can be expressed in terms of the basic parameters at the panel centers and 
network edges by equation (C7) yielding 

0 = Aj {SP ik } Nx 6 x k (D5) 


A similar relation holds for the contribution to 0 from the source distribution; but with 
the A^ for the source panels known from the boundary conditions of equation (D2), it 
contributes to the right hand side of the boundary condition equation (D3). The potential 
for a panel edge from the source distribution is given in equation (8.6). The contribution to 
the entire subpanel or panel is found by evaluating and summing equation (8.6) for each edge 
and at both endpoints proceding around the panel in a counter-clockwise direction in the 
same way as for the doublet panel. 

The pilot code computes the coefficients of the A^ in equation (D5) for each panel at 
all the control points before proceding to the next panel. A panel is selected and the control 
points for all the panels are transformed into the coordinate system of the subpanel by 
equation (Dl). Then a test is made to determine if the selected panel lies within or is cut 
by the Mach cone emanating upstream from the control point. The aerodynamic influence 
coefficients are not computed for those panels lying outside of the mach cone. The contri- 
bution from the selected panel is computed and stored for each control point. Then the 
program takes the next panel in order, repeats the procedures, and sums the results with 
the previous calculations for each control point. 
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Similar procedures are followed when velocity components are computed at the control 
points. The velocity components for the quadratic doublet distribution are given in equation 
(6.1) for the six parameter spline and in equation (7.4), (7.6), and (7.9) for the nine para- 
meter spline, in which the line vortices have been removed. The result of the computations 
of the boundary conditions in equation (D2) and (D3) lead to a set of simultaneous equations 
of the form 


N d N s 

S A ik*k = 2 B ik (U.fi k ) i = 1,2, . . . N d 
k= 1 k = 1 

to be solved for the doublet parameters X^. Here, and N g are the numbers of doublet 
parameter and source parameters, respectively. Once the X^ are known, a similar procedure 
is used to compute the velocity components and the pressure at each control point. 

For the combined source and doublet panels with the boundary conditions of equations 
(D2) and (D3), the tangential perturbation velocity components on the outer surface can be 
computed by differentiating the doublet strength on the panel. With the transformation 
from reference coordinate to panel coordinate system given by equation (Dl), the tengential 
velocity vector V r in the reference system is given in terms of the subpanel tangential velocity 
Vpby 


V r = KHA 1 } T lA 2 } T V p 


where the T denotes the transpose of the matrix. The perturbation velocity is used to compute 
the pressure coefficient using equation (3.10) or the isentropic relation. 

When the configuration also contains superinclined networks representing inlets, the 
procedure as outlined above is followed except the boundary conditions on the interior or 
downstream side of the network are 


0 = W . n = 0 

Since there are generally subinclined panels in the zone of dependence of the superinclined 
network, these boundary conditions yield relations similar to equation (D6) with the left 
hand side containing unknown source parameters as well as doublet. 

The velocity potential for superinclined doublet panels is given in equations (9.12) and (9.16) 
with the functions Qj and wq derived in Appendix A, section A9. The velocity potential for 
source panels is given in the equation following equation (9.16). The perturbation velocity 
components for the superinclined doublet panels are given in equations (9.30) and (9.31) and 
are used only in the nine parameter spline, since the line integrals of the doublet strength 
have been dropped. The velocity components from the superinclined source panel are pre- 
sented at the end of section 9.3. 
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APPENDIX E: 

RELATIONSHIP OF SUPERSONIC PANEL INFLUENCE COEFFICIENTS 
TO SUBSONIC PANEL INFLUENCE COEFFICIENTS 


The subsonic panel influence coefficients (ref. [3] ) were obtained by evaluating the 
integrals of the form 


H(M,N,K) 



(£-x)M-1 ( 77 — y )N 1 d%dri 
R k 


(El) 


and 


where 


F(M,N,K) = 


J 


(g-x)M-l (77-y )N~ 1 dg 
R k 


^ \/(£-x) 2 + Oj-y) 2 + h 2 


(E2) 


Here 2 is the panel surface (assumed to be flat for present purposes) and L is a typical edge. 
Also (x, y, h) are the coordinates of the field point in the local panel coordinate system and 
(f , 77 , 0) are the corresponding coordinates of the integration point. 

The integrals in equations (El) and (E2) can be formally defined for a subinclined panel 
in supersonic flow for which R becomes the hyperbolic distance defined by 


R ~\/(£-x) 2 - 0?-y) 2 - h 2 (E3) 


The H integrals satisfy three recursion relations obtained in the same manner as for subsonic 
flow. We have 


H(M, N, K-2) = H(M+2, N, K) - H(M, N+2, K) - h 2 H(M, N, K-2) (E4) 


4 

(K-2)H(M, N, K) = (M-2) H(M-2, N, K-2) - 2 vg F(M, N-l, K-2) (E5) 
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(E6) 


4 

(K-2) H (M, N, K) = - (N-2) H (M, N-2, K-2) + 2 vt F(M, N-l, K-2) 

1 

where the summation is over the panel edges. These recursions can be combined to produce 
the following, more useful recursions: 


4 

(M-N-K) H (M, N, K) = h 2 (M-2) H (M-2, N, K) - h 2 2 vt F(M-1 , N, K) 

1 

4 

+ 2 a F (M, N, K) 

1 


(E7) 


4 

(M-N-K) H (M, N, K) = -h 2 (N-2) H (M, N-2, K) + h 2 2 Vri F(M, N-l , K) 

1 

4 

+ 2 a F (M, N, K) 

1 


(E8) 


4 

(K-2) h 2 H (M, N, K) = (M+N-K-2) H (M, N, K-2) - 2 a F (M, N, K-2) (E9) 

1 


Here, a v £ and are defined as in the subsonic case, i.e., v = vr\) is the exterior normal to 

the edge K and a = v . (£ - x, r\ - y). 

The F integrals also satisfy three recursion relations obtained in the same manner as 
in subsonic flow. We have 


F(M + 2, N, K) - F(M, N + 2, K) - h 2 F (M, N, K) = F(M, N, K - 2) (El 0) 


vgF (M+l , N, K) + u v F (M, N+l , K) = a F (M, N, K) 


(Ell) 
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-(M-l)F(M-l, N, K-2) v n + (N-l) F (M, N-l, K-2) v% 

+ (K-2) F (M+l, N, K) Vfj + (K-2) F(M, N+l, K) ^ = E (M, N, K-2) (El 2) 

where 


E(M,N,K) = 


(g-x) M-1 (i?-y) N 1 
R k 


(E13) 


where the numbers 1 and 2 denote the endpoints of the subpanel edge. Because of the 
limited number of terms required for the flat panel supersonic case without line vortex terms 
no attempt has been made to recombine (E10), (El 1), and (E12) into more useful recursion 
relations. 

In order to evaluate the fundamental integrals H(1 ,1,3) and F( 1,1,1) it is useful to make 
the following definitions. Let 


Then 


and 


and 


b = - ^2, g2 = a 2 _ bh 2 , and £ = ^(ij-x) + 


R 2 = (g 2 - C 2 )/b 


F(l,l,l) =-g- 



dC_ 

R 


H(l, 1,3) =- 



,ad£_ 


R(a 2 -8 2 ) 


(El 4) 


(El 5) 


(E16) 


(E17) 


For linear source and quadratic doublet (without line vortex terms) influence coefficients 
F(l,2,l) and F(2,l ,1) are also required, and these may be obtained via equations (El 1) and 
(El 2). Upon considerable analysis we obtain the following evaluations: 
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hH (1,1,3) = 


0 , h = 0 (point not on panel) 

IT sign (hv£) , Rj = R 2 = 0 (wedge region) 

| tan - 1 £ haF j / (R j R 2 +h 2 F 2 )J , otherwise 


(El 8) 


where 


Fl 


( (*lR2- 
X (Ro 2 - 


-«2 Rl)/g 2 . h > 0 

(R 2 2-R 1 2)/( Cl R 2+ e 2 R 1 ) ; b < 0 


(E19) 


F 2 = 


(bRjR 2 + C}C 2 )/g 2 - b > 0 

(g2 _ j 2 _ C2 2 )/ (bR i R 2 - « 1 fi 2 ) b<0 


(E20) 


and the subscripts 1 and 2 on R and £ designate the quantities evaluated at the two ends of 
the side. Note that 

F 2 2 + b Fj 2 = 1 (E21) 

We also have 


F (1,1,1,)= 


Tr/V'b , Rj = R 2 = 0 
•— — tan~l (\/bFj/F 2 ) ; otherwise 

Vb 

sign (yjj) ^ / V^ibi" R 1 + I ® 1 1 


V^bl \/|bT R 2 + I £21. 


b > 0 


b < 0 


(E22) 


e [ 1 - b e 2 /3 + b 2 e 4 /5 - b 3 e 6 /7 + ...], F 2 > > V /1R I I 
where e=Fi/F 2 
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and 


I -[^(R 2 -Ri) +a^F (1,1,1)] /b ; b¥=0 

[-^(R 2 -Rl)Rl R2 + ^2 r 1 (^l-y)- c 1 r 2 0?2-y)] /g 2p 2 < E23 > 

- av n e 3 [l/3-be 2 /5 + b 2 e 4 /7 + ...] ; F 2 >> v /7bi |Fj| 


and 

F(2,l,l)= -vy (R 2 -R!) +a^F(l 1 l,l)-2^p T? F(l,2,l) (E24) 


The relationship between the calculations of equations (1) through (24) and those of 
Appendix A will now be established through the following identities: 


2 = v n (£ - x) + 1 >£ (t? - y) = - v v s m - vg, = v% (s m /m - s) 

dC = ^(l/m 2 - l) ds 
X = -v v lvj; 

m = ~v%/v v 

x m =-a/^ 

= (aoyj - /b 
t = (-aj'| + £i' TJ )/b 
Ql = -hH (1,1,3) 
w o = — ^ F (1,1,1) 

XWj + x m w 0 = ^ F (2,1,1) 

W] = ^ F (1,2,1) 

A 

x m = - a /^7? 
w D =v r\ F (1,1,1) 

w i ~~ v n F(i,2,l) 


(E25) 

(E26) 

(E27) 

(E28) 

(E29) 

(E30) 

(E31) 

(E32) 

(E33) 

(E34) 

(E35) 

(E36) 

(E37) 
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Figure 1. —Subpane f Structure of Basic Four Point Panel Showing The 
Five Planar Areas and Eight Triangular Subpanels 



155 


/ 

I 


x o*yo* z o 



Domain of dependence 

for x o» yo* z o 



of Xo, y 0 » *o 


Figure 4.— Domains of Influence and Dependence of The Point 
x O> V0' z o m Supersonic Flow 
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Figure 7. -Volume And Surfaces For The Uniqueness Proof From 

Boundary Conditions on Superinclined Surfaces in Supersonic Flow 



Figure 8.— Cross Section at y q = y of The Volumes and Surfaces For 
Deriving The Integral Equation For Subinclined Surfaces 
in Supersonic Flow 
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Figure 9.— Cross Section at y = 0 of The Volumes And Surfaces For 
The Uniqueness Proof For Subinclined Surfaces in Super- 
Supersonic Flow 



Figure 70.— Cross Section of Volume at y = Constant And Surface 
For Calculating Velocity Neglecting Edge Vortices 
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Figure 1 2.—Subpanel Structure of Basic Four Point Panel Showing 
The Five Planar Areas And Eight Triangular Subpanels 


161 


Source/Analysis 



Figure 13.— Schematic Location of Control Points And Corresponding 
Values of Source Singularity Strength For Analysis 
Networks 


Doublet/Analysis 



Figure 14.— Schematic Location of Control Points And Corresponding 
Values of Doublet Strength For Analysis Networks 

{Doublet/Design {Doublet/Design 




Figure 1 5.— Schematic Location of Control Points And Corresponding 
Values of Doublet Strength For Design Networks 





Figure 16.— Illustration of Mach Wave Propagation And Reflection of Interior Perturbations 



x »y (Projection of point x,y,z on panel) 


Figure 17.— Region of Integration For Subinclined Pane! Edge 







a. Mach cone intersection inside panel 




b. Mach cone intersecting one side of panel 



panel corner 
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d. Mach cone intersecting panel 
corner cone and opposite side 


Figure 20.— Examples of Mach Cone Intersections With Superinclined Panel Aligned 
Normal to The Free Stream 
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y e .z e 

Control point projection 


Figure 21.— Illustration of Panel Edge Coordinates With Pane i Edges Specified 
Counterclockwise on Exterior or Downstream Side of 
Superinclined Surface 



Figure 22.— Relation Between Panel And Edge Coordinate Systems of 
Superinclined Panel 
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Figure 23-Transformation From Reference Coordinates x r ,y r ,z r To 
Compressible Coordinates x c ,y c ,z c 



Figure 24.— First Angle of Rotation For Panel Coordinate System 



Figure 25.-lllustration of Second Angie of Rotation For Subsonic 

Flow or of Oblique Transformation For Supersonic Flow 
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Figure 26.— Comparison of Pressure Coefficient on Cones at Zero 
Angie of Attack With Exact Theory 



Figure 27.— Comparison of Three Methods for Computing Pressure Coefficients on 
10° And 15° Half -angle Cones at Zero Angie of Attack as a Function of 
Mach Number 
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Figure 28.— Comparison of Pressure Coefficients Computed by Zero Nc 
Zero Normal Velocity Boundary Conditions 
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Figure 29. -Comparison of Four Methods For Computing The Pressure Coefficient on The Spindle r = 0.08x2 (1-x) 







Azimuth angle 8 

Figure 32.— Comparison of The Pressure Distribution by The Present Method With The 
Van Dyke Second Order Theory For The Elliptic Cone Having a Maximum 
Half Angle of 14° And Fineness Ratio of 0.532 
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Van Dyke [21] second order theory 



Figure 33.— Comparison of The Pressure Distribution by The Present Method With The 
Van Dyke Second Order Theory For The Elliptic Cone Having a Maximum 
Half Angle of 75° And a Fineness Ratio of 0.3 
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Figure 34,-Comparison of The Pressure Distribution by The Present Method With The 
Van Dyke Second Order Theory For The Elliptic Cone Having a Maximum 
Half Angle of 18.4° And a Fineness Ratio of 0.33 
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Figure 35.— Comparison of The Pressure Distribution by The Present Method With The 
Van Dyke Second Order Theory For The Elliptic Cone Having a Maximum 
Half Angle of 30° And a Fineness Ratio of 0.2 



Figure 36.— Comparison of The Pressure Distribution by The Present Method With The 
Van Dyke Second Order Theory For The Elliptic Cone Having a Maximum 
Half Angie of 30° And Fineness Ratio of 0.2 With Refined Paneling And 
Modified Pressure Calculation For Near Stagnation Conditions 
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Figure 37.— Comparison of Pressure Distributions on an Elliptic Cone Resulting From 
Mass Flux and Velocity Type Boundary Conditions 
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Figure 38.— Comparison of Pressure Distribution From Pianar Doublet Pane / Method 
With Linearized Theory Fora Yawed Delta Wing 
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Figure 39,-Comparison of Pressure Distribution on Parabolic Cambered Wing From 
Planar Doublet Panel Method With Exact Linearized Theory 
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Figure 40.— Comparison of Downwash Computed From Planar Doublet Design Method 
With The Actual Wing Slopes on The Line N = 1 For Cambered Wing in 
Figure 39 
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Figure 42.— Comparison of Pressure Distribution on Symmetric Thick Wing Having Same 

Upper Surface Slopes as in Figure 39 With The Exact Linearized Theory Solution 
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Figure 43.— Illustrating Advantage of Surface Paneling With Exact Boundary Conditions Over 
Planar Paneling With Linearized Boundary Conditions 
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Figure 44,-Coefficient of Pressure on Outer Surface of Inlet Nacelle by Source Panel Method 
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Figure 45-Coefficient of Pressure on Outer Surface of Inlet Nacelle by Combined Doublet 
Source Pane! Method 




Figure 46. -Paneling on Carlson Wing 2T 
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Figure 47.— Comparison of Pressure Distribution Computed From Combined 
Doublet-source Panel Method With Experiment on Strips 2 and 4 
of Carlson Wing 2T 



Figure 48.— Comparison of Pressure Distribution Computed From Combined 
Doublet-Source Panel Method With Experiment on Strips 6 and 10 
of Carlson Wing 2T 
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Figure 49.— Paneling of Leading Edge Region of The Parabolic Arc Cambered Wing 
in Figure 39 Which Showed Instability in The Solution 
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Figure 50-Comparison of Pressure Distribution on Row 1 of Upper Surface of The Parabolic 
Cambered Wing Using Three Different Splines And The Paneling of Figure 49 
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□ Least square spline (6 coefficients/panel) 

O Corner point continuous spline (6 coefficients/panel) 
A Almost edge continuous spline ( 8 coefficients/ panel) 
() Edge continuous spline (9 coefficients/panel) 

— Theory 
Row 2 

30° sweep angle 

5 % parabolic arc camber profile 


Spanwise variable Y 


Figure 51.— Comparison of Pressure Distribution on Row 2 of Upper Surface of The Parabolic 
Cambered Wing Using Three Different Splines And The Paneling of Figure 49 
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Figure 52.— Comparison of Pressure Distribution on Row 3 of Upper Surface of The Parabolic 
Cambered Wing Using Three Different Splines And The Paneling in Figure 49 
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Figure 53.- Examples of Successful And Unsuccessful Paneling of The Parabolic Cambered Wing 
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Figure 56.— Comparison of The Pressure Coefficient on The B-1 Forebody Computed From The 
Doublet Gradient With Experiment And With The Finite Difference Solution 



Upper fuselage line 


O Experiment (Ref. [28]) 

• Source & doublet panels with 
six parameter spline 
— Supersonic finite difference 
method 


CP 


-.2 

0 


Cp 


.6 

0 

.2 



• 

V 

• s 

# 

mil 

.<■ 

9V9HU 

■ 


• 


t 

0i 


• '* 1 
•/ • 

l 

- • i 


T 7T 1 

A * 

^ Q 


— 

a 



• } 
1 1 






Lower fuselage line 


Upper fuselage line 




? 




t 

• 

• 



« 


Lower fuselage line 


0 5 10 

Body station, meters 


15 


Figure 57.— Comparison of The Pressure Coefficient on The B-1 Forebody Computed From The 
Aerodynamic Influence Coefficients Using The Six Parameter Spline With Experiment 
And With The Finite Difference Solution 
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Figure 58.— Comparison of The Pressure Coefficient on The B-1 Forebody Computed From The 

Aerodynamic Influence Coefficients Using The Nine Parameter Spline With Experiment 
And With The Finite Difference Solution 
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Figure 59-Comparison of Pressure Coefficient on The Carlson Wing Computed by Sources Alone 
W ith Experimental Measurements at Span Stations 2 and 4 
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Figure 60.— Comparison of Pressure Coefficient on The Carlson Wing Computed by Sources 
Alone With Experimental Measurements at Span Station 6 and 10 
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Figure 61 .—Comparison of Pressure Coefficient on The Carison Wing Computed by Combined 
Doublet-source Panels Using The Nine Parameter Spline With Experiment at Span 
Stations 2 and 4 
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Figure 62.— Comparison of Pressure Coefficient on The Carlson Wing Computed by Combined 
Doublet-source Panels Using The Nine Parameter Spline With Experiment at Span 
Stations 6 and 10 
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Figure 64. -Geometry of Stream Surface Generated at The Intersection of Conical Shock 
With Plane Parallel to The Shock Axis 
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Figure 65.— Pressure Distribution on Upper Surface of Asymmetric Configuration Bounded 
by Surface of Figure 64 and x,y Plane Using The Nine Parameter Spline With 
Combined Doublet-source Panels 
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Figure 66.— Pressure Distribution on Upper Surface of Symmetric Configuration Bounded 
by Surface of Figure 64 And Its Reflection in The x,y Plane Using The Nine 
Parameter Spline With Combined Doublet-source Panels 
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Figure 67-Pressure Distribution on Upper Surface of Asymmetric Configuration 

Bounded by Surface of Figure 64 And x,y Plane With Additional Doublets 
Distributed Over The Mean Surface With Zero Normal Mass Flux Boundary 

Conditions 
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Figure 68. - Pressure Distribution Near Leading And Trailing Edge of Delta Wing With 
Sharp Supersonic Leading Edge Using The Nine Parameter Doublet Spline 
With Combined Doublet-source Panels 
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Figure 69.— Pressure Distribution Near Trailing Edge of Delta Wing With Sharp Supersonic 
Leading Edge Using The Nine Parameter Spline With Combined Doublet-source 

Panels 
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Figure 71.— Pressure Distribution Near Leading Edge of Conical Flow Stream Surface Using 
The Nine Parameter Spline With Combined Doublet-source Panels 
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Figure 72 — 
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Figure 73.— Paneling on Wing/Body of NASA Memo 10-15-58L With Two Axes of Symmetry 
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Figure 77 .— Pressure Distribution on The Wing of The Arrow Wing/body of Figure 76 With 
Experimental Values at 9%, 20%, and 35% Span Locations and Computed 
Results at 12.8%, 22%, and 32% Span Locations 
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Figure 78.-Pressure Distributions on Wing of The Arrow Wing/body of Figure 76 With 
Experimental Values at 50%, 65%, And 80% Span Locations And Computed 
Results at 52%, 61.5%, and 80% Span Locations 
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Figure 79. -Paneling on The Arrow Wing/body With The Twisted Wing 


Combined source and 9 parameter 
doublet spline 

M = 1.7 
a * 4° 



Figure 80.— Pressure Distributions on The Twisted Arrow Wing/body Combination of 
Figure 79 at 9%, 20%, And 35% of Span Locations 
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Combined source and 9 parameter 
doublet spline 


50% span 65% span 80% span 



Figure 81. -Pressure Distributions on The Twisted Arrow Wing/body Combination of 
Figure 79 at 50%, 65%, and 80% Span Locations 



Figure 82.— Pressure Distributions on The Body of The Twisted Arrow W ing /body 
of Figure 79 



to 

to 

u> 


Mach 1.70 




0 * tt/24 


Mach number -\PT 

0 * 5 n /24 


.* 


5L 


Figure 84. -Comparison of Pressure Distributions Along Two Azimuth Positions, n/24 And 5ir/24 , 
on The Exterior Surface of The Nacelle For The Flows Calculated With And Without 
Superinclined Network in The Presence of a Perturbation From Constant Upstream 
Source Panels ■ 
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Figure 85. -Comparison of Pressure Distributions Along Two Azimuth Positions, 3ir/8 And 
13n/25, on The Exterior Surface of The Nacelle For The Flows Calculated With 
And Without Superinclined Network in The Presence of a Perturbation From 
Constant Upstream Panels 
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Figure 86.— Comparison of Exterior Pressure Distribution on The Exterior Surface of The Nacelle 
Containing an Interior Superincl ined Network With The Lighthill Solution 
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Figure 87.— Top And Side Views of Paneling on The LES 216 Supercruiser 
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Figure 88.-Perspective View of The Paneling on The LES 216 Supercruiser 
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Figure 91. -Evaluation of wq For One Edge Endpoint at Points on The Pane! Plane 


Mach wedge region 
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Figure 93.— Evaluation of Qf on The Pane / For Panel Subsonic Leading Edges 




Panel No. Other center panel values of source strength used in 

determining source strength in a given panel 

1 2 , 4,5 

2 1 , 3 , 4 , 5, 6 

5 1 , 2 , 3 , 4 , 6 , 7 , 8, 9 

Figure 94-Panels Used in Determining Source Strength From Weighted 
Least Square Fit 



Figure 95-Illustration of Subpanels Into Which Basic Four Point Non-planar 
Panel is Divided 
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Figure 96.— Illustration of Points at Which Doublet Strength Values Are Used 
to Define The Quadratic Distribution on Each Subpanel 



Figure 97.— Illustration of How Continuity of Doublet Strength is Obtained 
at Subpanel Edges 
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• Denotes point at which doublet strength is to be determined 
(X) Denotes points at which doublet values are used with w k = 1 - M ( t - R k • C) 

X Denotes points at which doublet values are used with w k = 10^ 


b 



Figure 99.— Illustration of Determining Doublet Strength at Center Point of 
Panel Edge and at Panel Comer by Weighted Least Square Fit 
With Neighboring Center Pane / Values For Interior Panels. 
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• Denotes point at which doublet strength is to be determined 
(x) Denotes points at which doublet values are used with weight 10^ 

X Denotes points at which doublet values are used with w^ = 1 - IVlCR^ C ) 
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Figure 101 -Illustration of Least Squares Fit Applied to Finding Corner And 
Midpoint Edge Values of The Doublet Strength For Panels Near 
Network Edges 
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